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(54) Title: TRAJECTORY ANALYSIS RADAR SYSTEM FOR ARTILLERY PIECE 



(57) Abstract 




A low cost easily deployed, yet highly accurate weapon mounted trajectory analysis radar system uses optimized Doppler 
radar signals to detect the actual trajectory of a projectile fired from an artillery piece such as a 155 mm Howitzer (10). The actual 
trajectory data is used to determine an atmospheric model that may be used to aim the weapon for future firings. The system in- 
cludes a Doppler radar system (16) having a turret mounted antenna (22) that tracks the projectile, a digital signal processor (30) 
using an FFT to convert radar pulses to trajectory data and a data processing system (32) that analyzes the trajectory data to de- 
velop a true atmospheric model. 
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TpA.TFgTQRY fiMAT^YSIS p&HAP SYSTKM FOR ARTTTT.'BRY PIECg 

gQpvriah-fc Authorization 
A portion of the disclosure of this patent 
document contains material which is subject to copyright 
protection. The copyright owner has no objection to the 
5 facsimile reproduction by anyone of the patent document or 
the patent disdpsure, as it appears in the patent and 
Trademark Office patent file or records, but otherwise 
reserves all copyright rigrfts whatsoever. 

Baelcarounci of the Tnvention 
3^0 AS weapons systems have developed, there have 

been substantial increases in both range and accuracy. 
However, systems that provide the greatest accuracy are 
extremely expensive, cumbersome to deploy and subject to 
failure because of their complexity. A need thus arises 
15 for a comparatively inexpensive weapons system that is 
reliable and accurate, yet lightweight and easy to use and 
deploy. 

The trajectory of motion and hence the impact 
point of a nonguided projectile fired from an artillery 

20 piece such as a 155mm Howitzer is determined by a 
relatively small number of parameters that are fairly well 
understood. These include the projectile shape, gun 
barrel characteristics, initial velocity of the projectile 
and atmospheric conditions. Typically all of these 

25 parameters will vary slowly over time, but if they can be 
determined accurately from one firing they will remain 
reasonably stable for a second firing shortly thereafter. 
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For example, the gun barrel typically contains 
helical grooves or rifling that causes the projectile to 
spin and attain greater stability and predictahxlxty as it 
exits that gun barrel. As the grooves wear, the impact on 
the projectile will change and the flight characteristics 
will oradually change with time. 

Similarly, the exact initial muzzle velocxty 
will Change gradually with barrel characteristics with 
temperature and with the powder charge that is ^--^ J^- 
powder charge will be quite consistent wxthm a 
manufacturing batch so that the initial muzzle velo"^-- 
for two consecutive firings will be nearly the same 
Atmospheric conditions including air density and wind 
velocity at different altitudes will tend to vary more 
rapidly than some of the other . parameters . However, even 
atmospheric conditions will usually remain reasonably 
stable for many minutes at a time. 

Radar systems have'^een developed to improve the 
accuracy of projectile firings. These systems tend to 
follow one of two configurations. One is a muzzle 
velocity radar system (MVR) that accurately determines the 
initial muzzle velocity. The other is a 
measurement radar system that tracks the tra3ectory of the 
projectile from firing to impact. 

The muzzle velocity radar system is compact, 
lightweight and relatively inexpensive. However, it only 
provides accurate information as to one of the many 
parameters that determine the final impact point, namely 
initial muzzle velocity. This information significantly 
improves the accuracy of a firing, but leaves many 
important parameters to be approximated by other means. 
Examples of this type of system can be found m U.S. 
patent 4,837,718 to Alon and U.S. patent 3,918,061 to 
Elgaard. 
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The trajectory measurement radar systems use 
tracking radar systems with multiple sensor points located 
some distance from a weapon. These systems are large, 
complex and difficult to properly deploy, especially under 
battle conditions where time may be critical. Because of 
the complexity of these systems and the distances over 
which they must be deployed, their reliability is 
questionable. The high power radar signals that track the 
small projectiles from a substantial distance are subject 
to detection and tracking by enemy forces. However, 
because these systems track the complete trajectory of a 
projectile , they can be used to compile extremely accxirate 
estimations of all of the parameters that affect the 
accuracy and final impact point of a projectile fired from 
15 a gun. 

U.S. patent 4,679,748 to Blomqvist et al. 
discloses a system that monitors the trajectory of a 
guidable projectile. The j^ojectile has flight control 
surfaces that are controlled in response to actual 
tracking information to guide the projectile to a desired 
impact point. The antenna is located off axis from the 
trajectory to provide the required trajectory information. 
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'-iin Y ""^ Tnvention 

A relatively inexpensive, lightweight trajectory 
analysis radar system for an artillery piece in accordance 
with the invention includes a gun mounted radar antenna 
5 tracking the trajectory of a projectile, a radar system 
sending radar pulses to and receiving radar reflections 
from the antenna, a signal processor analyzing the radar 
reflection signals to produce a representation of the 
projectile trajectory at least to the peak of the 
10 trajectory and a data processing system analyzing the 
trajectory data to determine initial muzzle velocity and 
atmospheric data. This data can then be output to a fire 
control computer system to enable the impact point for a 
next firing to be more accurately predicted. 

An advantageous trajectory analysis algorithm 
operates in a feedback loop using piecewise linear 
representation of an atmospheric model in selected 
elevation increments. The iodel is used to calculate a 
trajectory that is compared to the actual trajectory. 
Errors- bet*:*en the calculated and actual trajectory are 
used to update the atmospheric model for a next 

calculation iteration - 

Test simulation results suggest that two 
iterations will typically produce an atmospheric model 

25 that represents actual atmospheric conditions with 
sufficient accuracy. Atmospheric conditions may be 
assumed with adequate accuracy to be the same for both the 
rising and falling portion of the trajectory. Analysis of 
only the first half trajectory from firing to the peak is 

30 therefore all that is necessary to create an adequate 

atmospheric model. 

The predictable and relatively constant 
- projectile physical and flight -characteristics enable the 
radar system to be optimized with respect to size and 
35 power requirements, while minimizing the possibility of 
the tracking radar signals being detected by enemy radar 
units and used to determine the location of a weapon. 
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For example, mounting of the radar antenna on a 
recoilless portion of a gun turret with an approximate 
line of site along the gun barrel assures that the antenna 
will face the relatively large radar cross section of the 
trailing end of a projectile. The expected flight 
trajectory can be used to adjust radar power in accordance 
with increasing distance from the antenna and the phasing 
of transmitted pulses can be selected to prevent reflected 
pulses from interfering with transmitted pulses as the 
distance of the projectile from the antenna changes. Any 
required elevation control over the antenna can be limited 
and easily predicted in advance, while azimuth control may 
or may not be necessary, depending upon the requirements 
of any given application. The predictable velocity and 
15 position of the projectile can be used to eliminate the 
detection of false targets and thus improve the accuracy 
of the radar flight" trajectory data. 
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A better understanding of the invention may be 
had from a consideration of the following Detailed 
Description, taken in conjunction with the accompanying 

5 drawings in which: 

Fig. 1 is a perspective view of an artillery 
piece having a trajectory analysis radar system in 
accordance with the invention; 

Fig. 2 is a block diagram representation of a 
trajectory analysis radar system for a weapons system in 
accordance with the invention; 

Fig. 3 is a flow chart illustrating a computer 
program for analyzing trajectory data to derive 
atmospheric data in accordance with the invention; and 

Fig. 4 is a block diagram illustration of a 
method of analyzing projectile trajectory data in 
accordance with the invention. 
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p^va^ icH Deserip feion of t->te InvenHion 
Referring now to Fig. l, a weapons system 10 in 
accordance with the invention includes an artillery piece 
in the form of a 155mm Howitzer 12 having a track 
5 transport system 14 supporting a recoilless gun turret 16 
having a recoiling, rifled g-an barrel 18 mounted thereon. 
Mounted on the artillery piece 12 is a trajectory analysis 
radar system 20 having shown in Fig. 1 only an antenna 22. 

Antenna 22 is used to both transmit and receive 

10 radar frequency pulses for tracking a projectile that is 
fired from the barrel 18. Antenna 22 is preferably 
mounted on a recoilless portion of turret 16 and directed 
approximately along the line of sight of barrel 18. 
conventional azimuth and elevation motion controls are 

15 provided for antenna 22, which may be constructed to 
transmit a narrow .angle beam that tracks a reasonably 
predictable trajectory of a fired projectile. If the beam 
width is sufficient to assdre tracking of a projectile 
through all atmospheric conditions, it may be possible to 

20 eliminate azimuth motion control for antenna 22 in some 

circiimstanceis . 

AS a first shot is fired, the radar system 20 
tracks the projectile after it exits barrel 18 until it 
reaches the apex point of its trajectory. Tracking beyond 

25 the apex point is possible, but the information gathered 
tends to be redundant. At the same time, the radar power 
must be increased as the distance of the projectile from 
the antenna 22 increases. The increased operating time 
and the increased power would significantly increase the 

30 probability of enemy detection without significantly 
increasing the accuracy of future rounds. 

Enemy detection can be further minimized by 
having one weapon 10 fire a test round and then 
communicating the collected - meteorological data to other 

35 units via secure radio transmissions either directly or 
indirectly through a fire control center. 
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once data has been collected for one or two test 
rounds, further tracking would typically be discontinued 
until a significant change in the azimuth or peak 
elevation of the trajectory occurs or a significant period 
5 of time passes that might suggest a change in 
meteorological conditions along a contemplated trajectory. 
Even when full tracking is not used, the radar system 20 
would typically be operated under low power for a short 
time of typically less than one second to develop muzzle 
10 velocity information for each round. As each round is 
fired, the latest data for initial muzzle velocity for the 
same type of ammunition is used to predict nuzzle velocity 

for the next round. 

Referring now to Fig. 2, the trajectory analysis 

IS radar system 20 includes a digital signal processor 30 
which converts ra^dar data to a numeric coordinate data 
format based upon the center line of the transmitted radar 
signal, a' fire control dal^ processor system 32 and a 
doppler effect radar system 34. A main system data bus 40 

20 couples data processor system 32 to various functional 
units within the radar system 34 as well as to general 
fire control functions of weapons system 10 as represented 
by fire control operations 36. These fire control 
operations would typically include elevation and azimuth 

25 positioning of turret 16 and barrel 18 of weapons system 
10. Data processor system 32 is also coupled by a 
communication data bus 42 to digital signal processor 30 
and to a data output communication channel 44 which might, 
for example, connect to a secure digital radio 

30 transmitter /receiver. 

in preparation for firing an initial round, fire 
control data processor system 32 calculates an appropriate 
trajectory using the best available meteorological and 
initial muzzle velocity data and commands fire control 

35 operations 36 to position barrel 18 in accordance with the 
calculated data. The initial round can be either a first 
attempt at an effective round on target or a purely test 
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round that is fired at a high elevational angle in the 
direction of the target. A high elevation purely test 
round would have the advantage of providing meteorological 
data over a maximum altitude range and would be exploded 
shortly after reaching the apex of the trajectory. The 
highly accurate meteorological data derived from the test 
round could then be utilized to calculate with high 
accuracy the required trajectory for an effective round on 
target. 

in either event, the barrel 18 is positioned and 
a round is fired. An accelerometer 50 or similar device 
such as an acoustic wave detector detects the firing of 
the weapon system 10 and generates in response thereto a 
firing signal which is communicated to fire control data 
15 processor system 32 to establish time zero with respect to 
the firing of the. round. 

AS soon as the firing of a round is detected, 
data processor system 32 coifm'ands frequency synthesizer 52 
and transmitter unit 54 within radar unit 34 to become 
operational and begin transmitting radar pulse signals 
through a circulator 56 to antenna unit 22. Initially the 
projectile will be quite close to the antenna 22 and both 
the transmitted power and pulse width of the radar signals 
can be relatively low. As the distance of the projectile 
25 from the antenna 22 increases, the transmitted power and 
transmitted pulse width can both be increased to maintain 
adequate signal to noise ratios in the energy pulses 
reflected from the projectile back to the antenna for 
receipt thereby. In addition, as the projectile follows 
30 its normal trajectory, the phase and repetition rate of 
the transmitted radar pulses can be varied to assure that 
the reflected radar signals are received at a time window 
between transmitted radar pulses and to assure that a 
current pulse is not being transmitted at the same time- 
35 that a reflection from a previously transmitted pulse is 
being received. 



20 
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In a preferred embodiment, the radar pulses have 
a frequency in the 16-17 GigaHertz range with a pulse 
width of 0.5-2.0 microseconds. A pulse repetition 
frequency of 50-100 kHz results in a new pulse being 
5 produced every 10 to 20 microseconds. The initial power 
may be as low as l watt and is increased to as much as 2 00 
watts as the projectile reaches its apex. The antenna 
preferably has a gain of about 39 dB with an elevation 
beam width of approximately 1.70 degrees and an azimuth 
10 beam width of approximately 1.70 degrees. -It will be 
appreciated that other suitable frequencies and parameters 
could be used as well and can be optimized for any given 
situation. 

AS the weapon 10 is fired the firing of the 

15 charge creates an initial ionization zone around the 
barrel 18. Durii}^ this initial ionization period, the 
transmitter unit 54 goes through a warmup period and 
reaches the full commanded flower and the digital signal 
processor 30 ignores any return signals. After 

20 approximately 100 milliseconds this initial ionization 
dissipates and data processor system 3 2 issues a command 
over bus 42 to cause digital signal processor 30 to begin 
receiving and analyzing reflected radar signal data. 

The reflected radar signals received by antenna 

25 unit 22 include an azimuth difference signal, DIF A, an 
elevation difference signal, DIF E, and a sum signal. The 
DIF A and DIF E signals are communicated directly to an RF 
receiver 60, while the sum signal is communicated through 
circulator 50 to RF receiver 60. Receiver 60 mixes the 

30 three signals with a local oscillator signal from 
frequency synthesizer 52 with the resultant signals being 
communicated at an intermediate frequency of approximately 
500 kHz to an intermediate frequency receiver 62 for 
further amplification. IF receiver 62 mixes the signals 

35 with a coherent oscillator signal from frequency 
synthesizer 52 to obtain conventional velocity dependent 
Doppler signals and communicates the three signals to an 
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analog to digital converter at the input of digital signal 

processor 30. 

The azimuth and elevation difference signals 
provide indications of position error of the projectile 
5 from the center of line of sight of antenna unit 22. This 
information is in turn utilized during the course of 
flight of the projectile to command antenna position 
control 64 to reposition antenna 22 so as to maintain the 
projectile near the center of its line of sight. Antenna 
10 position control 64 may in turn provide back to data 
processing system 32 current actual antenna position data 
which can be combined with relative position data received 
from DSP 30 over bus 42 to permit calculation of the 
actual trajectory of the projectile relative to the outer 

15 tip of barrel 18. 

Digital signal processor 30 processes signal SUM 
with a 1000-2000 point fast Fourier transform (FFT) to 
obtain velocity data fromr'the doppler information of 
signal SUM. In addition, the time of occ\irrence of each 

20 reflection after transmission of the radar pulse signal 
corresponding thereto and the strength of the received 
reflected signal relative to the transmitted energy are 
utilized by digital signal processor 30 to generate range 
information that is converted to position information in 

25 an XYZ coordinate system utilizing the muzzle tip of 
barrel 18 as the origin. 

The three dimensional position and velocity data 
are calculated by digital signal processor 30 in response 
to the transmitted pulse repetition rate of approximately 

30 50-100 kHz. This high frequency data is smoothed by 
digital signal processor 3 0 with any data points being 
completed by interpolation of data on either side thereof 
and utilized to generate data points at a much lower 
frequency of approximately 100 Hz which are communicated 

35 to data processor system 32. Data processor 32 stores 
these data points for later trajectory analysis and also 
utilizes the data points on a real time basis to control 
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the positioning of ««nna 20 through 

control 64, to control th. power of th. tr«.«.itt.d radar 
i^^Ls, «.d to control th. pul.. repetition rate and 
pul.. duration of th. transmittal radar signals so as to 
optimize the efficiency of radar unit 34. As is 
eonv«,tion.l, digital signal processor 42 
data processor .yst«. 32 sampled data point 
a rate of approximately 100 ti»«i per second. The data 
Lr each sampled data point for the =-rv.d «a-="- 



XO includes three dimensional ^^^^^^ ^ ^^"^ 
dimensional position data, and a P^^^^^^.^^^J 
indicating the probability that the ^^'^^'^^^^J^''^'' 
point falls within a selected range of accuracy, 
information concerning the approximate ^^^^^^^ 
X5 position and velocity of the projectile can be utxlxzed by 
digital signal processor 30 to eliminate erroneous, noxse 
induced false detections of projectile posxtxon or 
velocity to improve the pro&bility and effectiveness of 
detecting true position and velocity. 
,0 If the data processor system 32 is fast enough 

to iceep up with communication and antenna positioning 
reguirements and still have additional processing txme, xt 
immediately begin analyzing the trajectory data 
received from digital signal processor 30 to determine 
25 initial muzzle velocity and accurate meteorological data. 
However, more typically, the data processor system 32 xs 
an Intel 80386 based microprocessor system that does not 
have sufficient speed and capacity to do both functions 
simultaneously. Furthermore, it is sufficient that the 
3 0 analysis data be available within a few tens of seconds 
after the projectile reaches its apex and there is no need 
for simultaneous computation of both antenna positioning 
-control and trajectory analysis. - 
The first 32 valid data points that are detected 
35 after firing are used to compute projectile muzzle 
velocity by performing a least squares error fit of the 32 
data points to a straight line. Since the initial 32 data 
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points will commence about 100 milliseconds after firing, 
the straight line is then used to extrapolate the velocity 
from the 32 data points back to firing time zero to 
establish the initial muzzle velocity. This value is 
5 stored and utilized to calculate the desired trajectory 
for the next round. Typically the muzzle velocity will be 
determined for each round and used to update information 
for the next round, even when a full half trajectory is 

not being tracked. 

The trajectory data, sampled at lO millisecond 
intervals, is then analyzed to establish the 
meteorological data including X and Y component wind 
velocity and air density at a plurality of different 
elevation points in the vicinity of weapon system 10. 
Meteorological data is typically calculated for sample 
points at 1000 foot elevation intervals. Polynomial 
approximation aAd moving average techniques are used to 
smooth the meteorologioari data resulting from the 
trajectory analysis and linear interpolation is utilized 
to estimate meteorological conditions between the 
meteorological sample data points which are stored for the 
different sample elevation points. In the coordinate 
system used herein X indicates a down range direction, Y 
indicates a cross range direction and Z indicates 
25 elevation. 

An algorithm for collecting the radar data and 
generating updated meteorological data is illustrated in 
accordance with the invention in Fig. 3, to which 
reference is now made. At steps 80 and 82 the actual 
projectile flight radar measurements are received and 
conventionally processed by digital signal processor 30 in 
real time during the actual flight of the projectile 



30 



At step 80 the digital signal processor 
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35 



analyzes the received radar data to produce sample points 
containing range, range rate, azimuth angle, elevation 
angle and probability of detection at a rate of 100 points 
per second. The data is determined relative to the center 
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line of the radar beam. Data processor system 32 keeps 
track of antenna position data received from position 
control 64 and later adds any offsets that result from 
motion of antenna 22. The data processor system 32 then 
continues to preprocess the converted data by selecting 
points having a probability of detection less than a 
certain threshold and replacing these points wxth 
estimated data derived from curve fitting a polynomial 
line through valid points on either side of the missxng 
points. The completed data is then prefiltered using 
integration and moving average techniques to smooth out 
any radar induced noise from the sample data- points. The 
three dimensional position and velocity information for 
the actual projectile trajectory is then used to derive 
accurate atmospheric information by in effect determining 
what wind and air density conditions would have caused the 

detected trajectory. 

once the projectile has reached the apex of its 
trajectory, the fire control data processor 32 is freed of 
its real time control functions and begins processing the 
stored actual trajectory data at step 84. The trajectory 
analysis proceeds at step 84 by establishing initial 
values for an atmospheric model or profile at selected 
small elevation intervals with respect to air density, p, 
cross range wind velocity v^, and down range wind velocity 
v.. Tne initial atmospheric profile is desirably made as 
accurate as possible and may be derived from any one of a 
variety of techniques. For example, conditions can be 
determined at ground zero and then initially it can be 
assumed that the ground zero conditions exist at each data 
point elevation level. Alternatively, the data derived 
from the last test firing can be utilized as the initial 
data. Another technique for initializing the atmospheric 
data may be to receive test fire determined data from a 
neighboring weapon system, from a fire control center or 
from other atmospheric data sources. 



SUBSTITUTE SHEET 



wo 92/19928 




PCr/US92/03377 



A closed loop iterative process to generate more 
accurate atmospheric data is then begun at step 8 6 by 
deriving a calculated trajectory for the fired projectile 
using the initial muzzle velocity determined from the 
5 radar measured actual trajectory of the projectile and the 
initial atmospheric profile Xq. The calculated trajectory 
produces position and velocity data at 0.01 second 
intervals corresponding to the 0.01 second data point 
intervals at which actual trajectory data is produced from 

10 the digital signal processor 3 0 processing of received 
radar data. At each of the corresponding data points, 
both position and velocity error for each of the three 
axial directions are calculated and this data is utilized 
to generate atmospheric error values for density, cross 

15 wind and down range wind at an altitude corresponding to 

the data point . 

Each of the six error parameters AX, AY, a2, 
AVx, AVy, and AV, is passed tl^ough a digital filter having 
the form 

20 (K/(l + ras))^a (1) 

where K is a gain parameter, r is a time constant, a is 
the current value of the parameter being filtered, Aa is 
the error value of the parameter being filtered and s is 
the LaPlace operator. In a preferred implementation of 
25 the invention, K and r have the values shown in TABLE I: 



K r 





0. 


05 


(in/s)/m 


0. 


1 


sec. 


AV. 


10. 


00 


(in/s)/(in/s) 


2. 


0 


sec. 


AY 


0. 


05 


(in/s)/m 


0. 


1 


sec. 


AVy 


10. 


00 


(m/s)/ (m/s) 


2. 


0 


sec. 


AZ 


0. 


01 


(kg/in')/ro 


0. 


1 


sec. 


AV, 


0. 


20 


(kg/m') / (m/s) 


1. 


0 


sec. 



TABLE I: Filter Parameter Values 
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Still Within step 86 the two filtered X 
direction error values for position and velocity, AXF and 
^VXF, are added to obtain a down range wind velocity error 
value AVWX. The two filtered Y direction error values for 
position and velocity, aYF and AVYF, are added to obtaxn 
a cross range wind velocity error value avwY and the two 
filtered vertical error values for position and velocity, 
tZF and AVZF, are added to obtain an air density error 

value Ap. *.v.^« 
These atmospheric model error values are then 

used at step 88 to update the current atmospheric model. 
During the first pass or iteration through loop 96, the 
current model is the starting model. The algorithm for 
updating the current atmospheric- model adds the error 
values to the values of the current atmospheric model. 
Typically the negative error feedback would be a 
subtraction, but it will be apparent that whether one adds 
or subtracts is merely a mister of the sign used for the 
error values. Those skilled in the art can use the proper 
combination of sign changes and addition or subtraction to 
implement a negative feedback loop and cause nhe 
atmospheric data model to converge toward the true 
atmospheric conditions. 

At tha end of a test firing sequence represented 
by get radar measurements 80 in Fig. 3. the data processor 
system 32 has received and stored data from DSP 30 at 10 
millisecond intervals defining the following parameters 
for the ballistic projectile: 

range 
range rate 
azimuth angle 
elevation angle 
. . probability of detection 

The data begins with the capture of valid data 
following the initial ionization interval and continues to 
at least the top of the ballistic trajectory curve. 



20 



25 



30 



35 
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The initial data format received from the DSP 3 0 
is then converted to three dimensional 6 degrees of 
freedom position and velocity format that places the 
origin at the projectile and produces for each 10 
5 millisecond sampled data point a set of position and 
velocity values relative to the firing point. The data 
produced is: 

position (down range) 
Py position (cross range) 
3^0 P. position (elevation) 

Vx velocity (down range) 
Vy velocity (cross range) 
V, velocity (elevation) 

This transformed data is computed from the 
15 original data for the corresponding sampled data points 
together with any offsets that occur as a result of 
positioning of antenna 22-- during the course of the 
trajectory. Relative antenna 22 position for each of the 
data points is received from antenna position control 64 
20 and stored by data processor system 32 in such a way that 
the antenna 22 position data can be correlated at this 
time with the radar analysis 10 millisecond sampled data 
point data. If antenna 22 is mounted on a turret of 
transport system 14, then any turret motion must be 
25 similarly considered. Using ordinary geometric 

relationships the position and velocity data can be 
calculated as follows: 



Px 
Py 



R . cos(E£) • cos(Az) (1) 
R . cos(E£) • sin(A2) (2) 



"30 p. » R • sin(E£) - ■ (3) 
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= R • COS(E£) • cos(Az) 

- R • E£ • sin(E£) • cos(A2) 

- R • A2 • cos(E£) • sin(Az) (4) 

Vy = R • cos(E£) • sin(A2) 

- R • E£ • sin(Et) • sin(A2) 

R • AZ • COS(Ee) • COS(AZ) (5) 



R • sin(E£) + R • E£ • cos(E£) (6) 



where R denotes range, R denotes range rate. El elevation 
angle, Et elevation angle rate, Az azimuth angle, and 
20 Az azimuth angle rate. The angular rates, Et and Az, are 
computued as the dif f eren9e- of a new measurement data 
point and the previous data point divided by the sampling 
time of 10 milliseconds. 

As the position and velocity data is calculated 
25 from the radar responses, a probability of detection value 
is determined from the signal to noise ratio of the 
reflected radar signal as calculated by isp 30. If the 
signal to noise ratio is less than or equal to 4 dB, the 
probability of detection is assigned a value of zero. For 
30 a s-n ratio greater than 4 dB and less than or equal to 
6 dB, the value is 0.1. If greater than 6 dB but less 
than or equal to 8 dB, the value is 0.3. If greater than 
8 dB but less than or equal to 10 dB, the value is 0.65. 
If greater than 10 dB but less than or equal to 12 dB, the 
value is 0.85. If greater than 12 dB but less than or 
equal to 14 dB, the value is 0.99. If the signal to noise 
ratio is greater than 14 dB~, a value'of 0.99"is assigned 
to the probability of detection. 
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As data for a sample point is processed, a- 
random numtoer is generated to create an event having a 
probability of occurrence equal to the probability of 
detection. For example, a random number between 0 and 1 
5 is generated and compared to the probability of detection • 
If the random number is less than or ec[ual to the 
probability of detection^ the event is deemed to have 
occurred and the corresponding sample data point is 
treated as valid. If the probability event does not 

10 occur, the sample data point is discarded and the missing 
point is determined from other points, either by a curve 
fit algorithm or by interpolation. A conventional least 
square curve fit of a polynomial to several points 
• surrounding the missing point is preferred. In this way 

15 data "it sampled points is used on a statistical basis in 
proportion to the probability that the data represents a 
valid trajectory point. 

An automatic weighing function is thus 
established where the data defining the sampled trajectory 

20 is weighed in favor of points having a higher probability 
of detection. However, points having a lower probability 
of detection are not completely ignored. 

Once three dimensional position and velocity 
sample data points have been calculated for each 

25 10 millisecond interval, the 6 sets of data are each 
independently smoothed by a suitable filtering technique 
such as use of integration or a moving average. The 
resulting smoothed data thus represents the actual, 
measured trajectory of the projectile at 10 millisecond 

3 0 intervals. This completes the processing of the actual 
trajectory data at step 82 in Fig. 3. 

At step 84 an initial estimate of the 
atmospheric model is established. The closer the initial 
estimate is to the true actual atmospheric conditions, the 

3 5 faster the determined atmospheric model will converge to 
the actual atmospheric conditions in response to 
processing using the actual trajectory data. One 
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technique is to initially assume that ^ 
exist Z all elevations of interest. -^^--^^/Z' 
initial estimate can be derived from a previous ^--J-J^^ 
the same or a neighboring location, from weather reports 
5 or forecasts, from observations or sxmply fro" ^ 
pre-established arbitrary estimate such as a typxcal axr 
density profile and zero wind velocity. 

The initial atmospheric data is established for 
air density; down range or X wind velocity VWX and cross 
10 range or Y wind velocity vwy". Temperature , T degrees, at 
" Te firing point is also established to -able future 
adjustment of air density and ballistic weapon parameters 
in response to temperature. However, the temperature xs 
assumed to be correct and is not modified in response to 
15 the actual trajectory data. The atmospheric --^-^^ 

stored in a table in small elevation increments begxnnxng 
with the firing -point as elevation, Z, equals zero. 
Typically data is stored at^n elevation correspondxng to 
e^h sample data point along the projectile trajectory and 
20 interpolation is used to find atmospheric data at 
elevations between these sample data points. Data could 
alternatively be stored for fixed elevation increments, 
increments of no more than 1000 feet are preferred. 

Linear interpolation is used to obtain 
25 atmospheric parameters at elevations between the stored 
elevation points. A more sophisticated interpolation or 
curve fitting could be used if desired, but linear 
interpolation is convenient and is sufficiently accurate 
to produce good results. 
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In the present model the atmospheric conditions 
vary only with altitude and for a given altitude are 
assumed constant and uniform throughout the range of the 

weapons system 10. 
5 The method of starting with an initial estimate 

of an atmospheric model and causing the model to converge 
toward an accurate true atmospheric model is analogous to 
a digital position feedback servo system as illustrated in 
Fig. 4. The converging process can be more clearly 
10 understood by examining the same process from two 
different perspectives as illustrated in Figs. 3 and 4. 

Pig. 4 represents a computer simulation of a 
method in accordance with the invention of developing an 
atmospheric model in the vicinity of a weapons system 10. 
15 At block 102 a current, accurate true atmospheric model of 
temperature, air. density, down range wind velocity and 
cross range wind velocity is assximed. 

At step 104 the projectile equations of motion 
are used to calculate a simulated actual trajectory that 
20 would occur given the true atmospheric model assumed at 
step 102. A 6 degree of freedom (DOF) mathematical model 
has been developed for all major weapons systems including 
the 155mm Howitzer weapons system which is representative 
of the weapons system 10. This model permits the 
25 assumption of any initial condition such as temperature, 
muzzle velocity, azimuth angle, elevation angle, 
projectile spin rate for each different aerodynamic 
projectile shape or other physical characteristics. The 
model is conventionally applied to an atmospheric model to 
30 calculate a projectile trajectory. One well known model 
is the BRL aerodynamic model from the Ballistic Research 
Laboratory at the Aberdeen Proving Grounds. 
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In any event, at step 104 the selected 
projectile model is applied to the atmospheric model 
established at step 102 to generate simulated sample data 
points at lO millisecond intervals representing in three 
S dimensions the position and velocity of a fired projectile 
such as a 155mm projectile. 

subsequently at step 106 the position velocity 
data of step 102, 104 is converted to a radar 
representation format and a radar model is applied to the 
10 calculated trajectory to superimpose simulated radar noise 
and probability of detection information. Some of the 
probability information is below the detection threshold, 
necessitating the estimation of missing points from 
surrounding data, just as would occur in the generation of 
15 true radar analysis data. Steps 102, 104 and 106 of 
Fig. 4 thus correspond to step 80 in Fig. 3. 

At steps 108 and 110 the radar simulated 
trajectory data is converted to a position, velocity 
format and smoothed by a prefilter step 110. The 
20 resulting data is a set of simulated dated points at 
10 millisecond intervals that simulates the output from 
process measurements step 82 of Fig. 3. 

An initial estimate of an atmospheric model is 
then established at step 84 in Fig. 3 and at a 
25 corresponding step 112 in Fig. 4. These initial values 
are used to initialize a current model data set at step 
114 of Fig. 4. The initialized current atmospheric model 
from step 84 of Fig. 3 is applied to a loop 96 and 
similarly the current atmospheric model from step 114 of 
30 Fig. 4 is applied to a loop 116 that corresponds to loop 

96 of Fig. 3. 

The loops 96 and 116 are functionally identical 
and operate on the current atmospheric model to cause the 
values represented thereby to converge toward the true 
35 correct atmospheric model represented by true atmospheric 
conditions in Fig. 3 and by atmospheric conditions 
simulated at step 102 of Fig. 4. 
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At Step 8 6 of Fig. 3 an aerodynamic model of the 
projectile is applied to the current estimate of the 
atmospheric model to calculate a ballistic trajectory. 
The calculated ballistic trajectory is then compared with 
5 the true trajectory data developed at step 82 and the 
position and velocity difference or error information can 
be used to calculate error values for the current 
atmospheric model. The error data added to the current 
atmospheric model values at step 88 to form a new, rough 

10 atmospheric model. The new, rough atmospheric model is 
smoothed ar step 90 and transferred to form a new current 
atmospheric model at step 92. 

At step 94 of Fig. 3, an exit criterion is 
tested. In the present instance it is determined that 

15 loop 96 has experienced only one iteration and control 
returns to step 86 to execute a new iteration with the 
smoothed new atmospheric model from steps 90 and 92 being 
used as the current atmosph^ic model. 

When the exit criteria is met at step 94, the 

20 loop 96 is exited and the then existing current 
atmospheric model becomes the assumed actual atmospheric 
model to be used in aiming calculations for subsequent 
firings of weapons system 10 and neighboring weapons 
systems • 

25 The loop 116 of Fig. 4 is functionally identical 

to loop 96 of Fig. 3 and begins at step 120 with the 
application of a projectile model to the current 
atmospheric model from step 114. Step 120 is identical 
to step 104 except that a different atmospheric model is 

30 used to calculate the trajectory. The trajectory from 
step 120 is subtracted from the calculated true trajectory 
at subtraction step 122 to produce trajectory error values 
that are applied to a transformation and smoothing step 
124 that generates atmospheric model error values aVWX, 

35 AVWY and Lp . 
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The transforming and smoothing function 124 uses 
the equivalent of position plus velocity feedback to 
filter and sum the trajectory error values to produce the 
atmospheric model error values. Although not explicitly 
5 shown, it will be understood that the atmospheric error 
values are determined for the elevations computed by block 
120 at 10 millisecond intervals. 

Transformation step 124 applies to each of the 
6 trajectory error positions and velocity error values a 
10 low pass filter function 130-135 represented in LaPlace 
notation by 

E* » (K*/(l + T*S))A* (7) 

where E* is one of the 6 trajectory error values, K* is a 
scaling or weighting value, r is a time constant and s is 

15 the LaPlace operator. 

An adder" 136 sums the weighted and low pass 
filtered X position and vel^ity trajectory error values 
to produce an X component wind velocity error value AVWX. 
An adder 138 sums the weighted and low pass filtered Y 

20 position and velocity trajectory error values to produce 
a y component wind velocity error value aVWY. An adder 
140 sums the weighted and low pass filtered Z position and 
velocity trajectory error value to produce an air density 
error value A/>. subtracter 122 and transformer 124 thus 

25 correspond to step 86 of Fig. 3. 

These filter equations are set forth below using 
the weighting values and time constants set forth in 
TABLE X. 

AVWX = AXF + AVXF (8) 
30 AVWY = AYF + AVYF (9) 

Ap = AZF + AVZF (10) 
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AXF « (0.05/(1 + O.lx s))AX (11) 

AVXF = (10/(1 + 2Vx s))AV, (12) 

AYF = (0.05/(1 + O.ly S))Ay (13) 

5 AVYF « (10/(1 + 2V^ s))AVy (14) 

AZF « (0.01/(1 + O.lZ S))A2 (IS) 

AV2F (0.2/(1 + IV, S))AV, (16) 



An adder step 142 corresponds to step 88 of 
Fig. 3 to add the atmospheric error data to the current 

10 atmospheric data to produce new atmospheric data that is 
smoothed at step 144 aria applied to the current 
atmospheric model at step 114 to replace the original 
current data with the new, smoothed atmospheric data. An 
exit substep 150 is then performed within step 114 to test 

15 for exit conditions. Upon exiting^ the current 
atmospheric data of step 114 becomes the data used for 
aiming calculations for weapons system 10. 

It will be noted that the adder 142 of Fig. 4 
corresponds to the update step 88 of Fig. 3, that 

20 smoothing step 144 of Fig. 4 corresponds to smoothing step 
90 of Fig. 3 and that transfer of the new atmospheric 
model to the current atmospheric model at step 114 of 
Fig. 4 corresponds to step 92 of Fig. 3. Test 150 of 
Fig. 4 corresponds to exit test 94 of Fig. 3. 

25 In computer simulations of the present 

invention, it was found that two iterations of loops 96, 
116 were sufficient to cause the current atmospheric model 
to adequately represent the true atmospheric model. A 
greater or lesser number of iterations could of course be 

30 used as the exit criteria. Alternatively, exiting could 
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be based on the magnitudes of the trajectory or 
atmospheric error values; i.e., when these values drop 
below a selected threshold. 

A trajectory analysis radar system siaulatxon 
5 program written in the Fortran programming language using 
the ACSL (Advanced Continuous Simulation Language) to 
implement the projectile equations of motion 120 and the 
low pass filter function 124 is set forth in APPENDIX A 
hereto. 

3^0 While particular arrangements of a trajectory 

analysis radar system have been shown and described for 
the purpose of enabling a person of ordinary skill in the 
art to make and use the invention, it will be appreciated 
that the invention is not limited thereto. Accordingly, 
15 any modifications, variation or equivalent arrangements 
within the scope, of the attached claims should be 
considered to be 'within the scope of the invention. 
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APPENDIX A 



1 



PniMTFD FROH NODE - VAX2 

AT TIME = .16-OCT-1990 15s07j03.ei 

riLF CREATION DATE - 16-OCT-1990 .15:05:26.96 

LAST REVISION DATE - 16-OCT-1990 15:06:09.19 

PLUCK SIZE OF FILE - 341 

COPYRIGHT (C) LEAR ASTRONICS CORPORATION 1991 
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Fll«« RADAR_COM18.rOR 



C IN'-LUDE FILE FOR RADAR COHPUTATIONS 
C names RADAR_COM18 . FOR 



O2-MAY-1990 



IMPLICIT NONE 

INTEGER ISEED/2234567/ 

C Radar parameter definitions and units 

! The dB unit of X is defined as 10*LOG10(X) 
! variables set by user ^-ia routine RAD^f^^ 



! RADAR POWER 
! ANTENNA GAIN 
I WAVE LENGTH 
! RECEIVER FREQ 
I NOISE FIGURE 
! TOTA1~LOSSES 
1 rrT t.ENGTH 
! WEATHER 



radar uower 

antenna gain 

wave length 

receiver frequency 

system noise figure 

sYSt9H» losses 

f ft length used 

clear (0) or 4 wm/Hr rain 



(1) 



variables coiiting from ^I^^^^^f^l^P^f--— 
PROJ^C^IAMETER s projectile diameter 
variables computed in routine RADARSET 



SECT TON 

EI.EV BEAKWIOTH 

AZIM^BEAMWIDTK 

RANGE RES 

RRATE'RES 

RCSMAX 

RCSMAXDB 

STN 0 

RCS^NORM 

STNDB 

PROB DETECTION 
ATTEHfJATION 



prc» iectile section area 

elevation beamwidth 

azimuth beamwidrh 

range resolution 

ranae rate resolution ^ 

maximum radar cross sectxon 



free space S/N at 1 NM 
normalized RCS ([0.1]) 
S/N ratio 

probability of detection 
attenuation due to weather 



for max RCS 



Units : 

W 

dB 

m 

Hz 

dB 

dB 

n/a 

0/1 



m**2 

rad 

rad 

m 

m/s 
m**2 

dBsm 

dB 

n/a 

dB 

n/a 

db/km 
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C 'Declare radar parameter types 

^^^^ onwFR ANTENNA GAIN , WAVT LENGTH 

55SdB^ : l5oB DETECTION : ATTENUATION 

X BAND(2) KU BSNDi2) , 5-|^S 1 

KX BANbfi) : Q_BAND<i) , VlBAND(2) 

W_BAWD(2) 

C Conanon block for radar parameters 

COMMON /RAP^J_i2"powFR ANTENNA GAIN , WAVE LENGTH 

SSSf''^ : IrOB^ETECTION , ISEED 

ATTENUATION . ^ BAHD , q^bInD 

K BAND , S^wfi^^ ' Q.BAND 

; , V-BAND r W_BAND 

Q Default values for radar parameters 

^^"^^ RADAR POWER , S^|?»*|„^^^So ' SoilB^'lliSSE 

I?|rBEAMWIDTH . WEATHER , SECTION 

*' 155:e.3 : 0.5|6 . 2 

:: oioaoi : o. - o.oies? 

X BAND / 8.E9 , 12. E| / 

KU BAND / 12. E9 , 18. E| / 

K BAND / 18. E9 , 27. E9 / 

KX BAND / 27. E9 , 40. E9 / 

Q BAND / 36. E9 , 5f / 

V-BAND / 46. B9 , 56. E9 / 

W^BAND / 56. E9 , 100. E9 / 

C Constants 
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r — — — 

r —— — — — — • 

c 

prAi. 



PI 

DEG2EAO 



convention 
X (deg) 



(nm) 
X (deq) 
y irad) 



, TWOPl _ 
, RAD2DEG 



DATA 



NM2K 

DEG2KAD 

RAD2DEG 



r TWOPI 

ftlHiiiiJ., ; Sir 



NM2M 

LIGHTSPEED 



NM2M , 
LIGHTSPEED / 

1853. . 
3.E8 / 



WLXmt TAR.COK18.FOR 



e VARIABLE PECLARXTION AH^ 



COMMON BLOCK FOR TAR.MAIN PROGRAM 



CHARACTER* 10 
INTEGER 

Ir 

1, 
REAL 

ir 
Ir 
Ir 
Ir 
Ir 
Ir 
Ir 
Ir 
Ir 
Ir 
1 r 
Xr 



ICOFILE, ICFILEr TMPFILE 



OPTION. 



RHOTAB 



IHAX, 
ISEED 



JHAXr 



RNSEEp/223457/ 

NCOST(40) 



(40) , VWXTAB 
RHOTABOMO) 
PHOUL ' ' 
RKOLL 
RHOPCT 
t'OSTVX 
LAYER, 
RHOSTEP 



(40) 
VWYTAB0(40) 



40) 
(40) 
(40) 
(40) 

VWXSTEP 



(40) r VWYTAB 
VWXTXB0(40) 
VWXOL (40) 
VWXLL (40) 

vwxpcrr (40) 

L5?!^2i'cisTV2MIN 
VWYSTEPr XRHO 



VWYUL 
VWYLL 
VWYPCT 
COSTVZ 



(12> 

(40) 
(40) 
(40) 

COSTVZTH 
RHOMIN 



VWXIlO , 
VWYHO. 
KXER, 
K-VER f 
rtr.ERr 



VWXAO 
VWYAO 
KVXERr 
KVYERr 
KVZERr 



VWXHlr VWXAlr 

VWYHlr VWYAl 
KAXER 

kayer 

KAZERr KRGERr 



COSTVZTHSL 



KRRER 
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COMMON / TARCOMl / 
1 

1, 
1, 
1, 
1, 
1, 
1, 
)., 
1, 
1, 

Ir 
1, 
1, 
1, 
1. 
1, 
1, 

DATA OPTION, RADFLG / 2 . 2 / 

DATA ICOFILE. ICFILE , TMPriLE /.'ICO', 

DATA LAYER / 1000. / 



OPTION, RUNFLG, RADFLG, RNSEED „^ 
SLAYER. NLAYER, MLAYER, ILAYER, KLAYER 
IMAX, JMAX, KMAX'^NCOST 
RHOTAB , VWXTAB , ^^J^l^ 

MHOTABO , VWXTABO , y"^TABO 

m\oVL , VWXUL , VWYUL 

tJHOLL , vwxLL , yyj^^j;- 

RHOPCT , VWXPCT . VWYPCT 

COSTVX , COSTVY , COSTVZ 

^•;i5i?EF^'^xl¥lp,''^STEP, RHOMIN 

;;^S§: ^^tS: ^iii: ^Si 

KXER, KVXER, KAXER 
KYER, KVYER, KAYER ,»bpd 
K7.ZR, KVZER, KA2ER, KRGER, KRRER 

ISEED 



'IC, 'TMP' / 



DATA (PHOTABO(I) .I;1'20) 

1/ i.22<5, 1.1095, l-221i' 
1, 0.7357, 0.6559, 0.5872, 
1 0.4112 0.3616, 0.3118, 
1, 0.1929, 0.1640. 0.1423 
1/ 

DATA VWXIlO, VWXAO, VWXHl, VWXAl 
DATA VWYHO, VWYAO, VWYH.U VWYAl 
DATA LVX , LVY , LV2 / 1 • , 1 • , 1 • 
DATA IMAX, JMAX, KMAX / 50, 20, 



0.9082, 
0.5218, 
0.2629, 
0.1206, 



/ 
/ 
/ 

20 



2200. 
2300 



.8184 
.4656 
.2259 
.1029 



10 
5 



5100 
3800 



20 
15 



DATA (RHOUL(I),I-1'20) 
1/ 1.3000, 1.2000, 
1. 0.8000, 0.75OO, 
1, O.8O00, 0.80OO, 
I. 0.8000, 0.80OO, 
1/ 



1.1000, 1.0000, 0.9000 

0.7000, 0.6000, 0.8000 

0.8000, 0.8000, 0.8000 

O.IOOOI 0.8000, 0.8000 



DATA (RHOLL^IJ,I-1,20^ 0.9000, 0.8000, 0.7000 
1 06000 0 5500 y.5000, 0.4000, 0.4000 
l' O'.IOOO: oioOO, 6.2000, 0.2000. 0.2000 
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I 0.0500. O.OSOO. 0.0500, 0.0500, 0.0500 

1/ 

DATA {VWXUL(I).I-1'20) ' 20*25. / 
DATA (VWXLL(I),I»1»20) ' 20* 0. / 
DATA (VWYUI,(I),I-lr20) • 20*20. / 
DATA (VWYI.I,(I).1-1.20) ' 20* 0. / 
DATA COSTV2THSL / 0.1 / 

DATA KXER, KVXER, K^xp . 0 05, Jg;; g; 

DATA KSERr K^ZERr KA2ER O.uir ^ 
DATA KRGER. KRRER / O-r 



PROGRAM TAR_MAIN18 



I 02-MAy-1990 



MAIN FRHGRAM FOR TRAJECTORY ANALYSIS RADAR (TAR) 

OPTION 



5 RUNFLG 



r 
r 
C 

c 

r 

r 
C 

c 
c 
c 

r 
r 



« 3 

« 1 

* 2 

« 3 

» 4 



PUN NOMINAL^AND riNAL{s) OHLY 
RUN PARAMETER GRID IN ONE LAYER 



RADFLG 



0 or (1) 
2 



mm ctosBD-LOOP estimation 

NO RADAR R>NI>9M 
RANDOM NOISE ACTIVATED 



,^ .^^r, CTM r«!T.' acSL DYNAMIC SIMULATION PROGRAMS 
THIS PROGRAM CONTROLS 'TAR.SIM.CSL AC5I* uxi««i 

- SPECIFY A ;«COM-ICriL| JO INITIALISE ACSL^RUN^^^^ ^ 

- DETERMINE RUN PARAMETERS (EG, wi^ SUBROUTINE 

- SET RON PAWmETERS THRU SET Ata^ *-! PARAMETERS 
: ^?Rli^ ^^rROr^SL^THRU^ETlACSL.PARAM SUBROUTINE 
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PARAKETERS PASSED TO TAF^SIM BEFORE EACH RUN ( SET^ACSL.PARAM) 



RUNFLG 
RADFLG . 
NLAYER 
LAYER 
RHOTAB 
VWXTAB 
VWYTAB 
LVX, LW 



RNSEED 



APFECT TERKT AND DATA LOGGING IN ACSL RUN 
RANDOM NOISE ACTIVATION AND SEED 
umliR Sf ATMOSPHERIC LAYERS TO BE RUN 

Wf^IGHTS OF VELOCITY ERRORS IN COST FUNCTIONS 



PARAMETERS RETRIEVED FROM TAR.SIM AFTER EACH RUN ( GET^ACSL^PARAM ) : 



Mr^AVFR 
CUSTVX 
C05TVY 
C05TV2 
NCOST 



l.AYHtn INDEX OF ( NOMINAL RUN ) APEX 

SSI g? 155?!! ^ S 1 ' 15 ffii 



C KEY ESTIMATION PARAMETERS DISPLAYED: 



r 

C 

c 

r 

C 
C 
C 

c 
c 
c 

C 

c 
c 

r 

c 

r 
r 
C 
C 
C 



RHOPCT 
VWXPCT 
VWYFCT 
D I MP AC 



rrRCENT ERROR IN RHO 
reRCENT ERROR IN VWX 

SSSilM^SNiE SSOR AT GROUND IMPACT 



CONSTRAINTS ON PARAMETER SPACE ARE: 



RHOUL, RTIOLL 

VWXUL, VWXLL 

VWYUL, VWYLL 

NOMINAL WIND BREAKPOINTS COMPUTED FROM A SMOOTH CURVE: 

VWXHO, VWXAO WIND X: AO IS CST WIND AMPLITUDE BELOW ALT HO, 

VWXHii VWXAl 

VWYHO, VWYAO 

VWYHX, VMYAi 



WIND Ys AO IS CST WIND AMPLITUDE BELOW ALT HO, 



C IMPORTANT ATMOSPHERIC LAYER INDEXES AND TERMINOLOGIES: 



C VARIABLE 

C 

C 

C LAYER 



SOURCE 



TAR MAIN 



DEFINITION 
SIZE OF EACH LAYER IN METERS 
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^^^Ul ^^iS^IH i^HllSI ffii"M"!l^VER, FOR SIM 

^ l^^ILE TArImAIN "COM FILE TO |TART | lAy|R igj^YER" 

^ MlS^i' Arg^^S-^^'ACSL liiiilB^^i ^?Ml?^i- TO RUN i "SLAVER- 

IMPLICTT NONE ^ . 
INCLUDE -TAR^COMXB.rOR- 

cllim inltiaiire ACSL 

RUNFLG^=^1^ I Always run Nominal first 

rIIIIII^Use"clefined run P*""*®^®ff - - — 

XOO CALL GET_USER_PAPJ\M 

g::::::;;;;;;^;^^^^ 

IF (RUNFLG.EQ.l) THEN 
EHO IF 

GOTO (lOOd, 2000, 3000) OFTIOH 



1000 
502 



PRINT 
FORHAT( 



(ix!70( J ,/rlX, 'Estimation Run ' r / r IX, 70 ( ' - ' ) ) 



SUBSTITUTE SHEET 



wo 92/19928 



PCT/US92/03377 



35 

cII"II''"Run*3-axis"cIose*" 

IMAX » 2 

Number of iterati'-^ts in DO loop 

rPINT * , ' Change number of iterations IMAX 7 IMAX 
FEAD(5,*) IMAX 
r — . Guess of rho 

REALM 5r*) (RHOTAB(I) .20) 
i; Guess of wind 

rnTMT * ' rnter fir«^- aness of VWXTAB and VWTTAB ?' ^ \, ^ ■ 

rlun lo: '^X?AB^ -?'IvwXTAB (I). 1-1,20) I Current Estxmate 

<-;;S$?jtB^'^'''^(^TAB (I). 1.1,20) I current Estimate 

READ( 5 , * ) ( VWYTAB( I ) , I«i . 20 ) 
10 rORMAT(lX,All,3X,JL0(F5.2,lX),:,2{/l5X,10(r5.2,lX))) 

f; Perform lat closed-loop run 

I « 1 

^ PRINT SXl^^^^^j^ .^Q^,^. J 'Iteration ',!) 

CALL RUN^ACSL(O) » l^oad Param/Run ACSL/Get Cost Function 

C Perform remaining runs 

DO I = 2, IMAX * Maximum number of iterations 

?^rRu"Ack(l) I l-oad Param/Run ACSL/Get Cost Function 

END DO 

C — — — Run Final(s) -♦•••♦♦•••♦•♦••^ 
iSS° FOwSt( ix! ? 0 (•-•.,/. IX , ' Final Run ' , / , IX , 7 0 ( ' - ' ) ) 
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RUNFLG = 3 

r 2050 CALL SET METS 
" CALL RUN_ACSLCO) 

. Select METS again 

C GOTO 2050 

PRINT • TAR> Rerun nominal ? {l«ye». ©.no) . 
READ(5,*1 I 

ihTiih^l''. ICOFTM) . r..«r. .t.«»P ..I— 

CALl. 2ZSVRS2(1, TMFFILe) 
END IF 

CALL GBT.ACSL PARAM ( ^ ^ , COSTW , COSTVZ , HCOST 

1; :^o?AB :^A8 i 

GOTO 100 



i:^:::;5;;:^«;:UKlsjis.i^!njs!..:.ii;.^:.5»^^^ 

RUNFLG = 2 

r First, Re-generate TMPFILB up to SIAYER 

' PRINT * , • TAR> Which layer to run grid ? 

»oT«T Bl ' Change VWKTAB - • , (VHXTAB(I) ,1-l.IIAYER) 
^orm. 81 • Change VWYTAB - • , (VWVTA8(1) ,I-1»IIAYER) 



SUBSTITUTE SHEET 



wo 92/19928 



PCT/US92/03377 



37 

Bl FORMAT(/A18,<ILAyER>ri0.3, ' ?') 

.''•iirRUNj^Ml) GenerSf TMPFILB 

C Grid run at layer '=5L.ArER 

SLAYEP = ILAYER 
KLAYER = ILAYER * 1 
NLAYEP • .1 

PRINT 91. IMAX, JHAX, KMAX 

Format?/ lx!^Ti^'''^.t.r IMAX, JKAX. KMAX for Grid Run .',315) 

PRINT 92, RHOSTEPr VWXSTEP, VWYSTEP 

i!^^'''*^ORSX$?7!x:'T)[$n^teriiSi?EP, VWXSTEP, VWYSTEP :-,3ri0.3, 
PRINT 03 RHOI,l.(KLAYEP> - VWXLL ( KLAYER ) , 5 55^^12 1 ' "^^^'^ 

DO 3200 1=0, IMAX-1 
DO 3200 J » 0, JMAX-1 
DO 3200 K « Or KMAX-1 

rr. /Twav r*x 1^ RHOTAB f KLAYER ) ^ RHOLL( KLAYER ) + RHOSTEP*I 

3 200 CONTINUE 
GOTO 100 ' Select OPTION again 

END 

SUBROUTINE RUN_ACSL(HEWRUN) 

i SUBROUTINE TO INTERFACE WtTH ACSL DYNAMIC SIMUIATION TAR.SIM PROGRAM: 

r PFRFORM AN •ACSI.> RESTOnf ' FILE (STARTING POINT) 
C : lIaE RELEVANT OPTIMIZATION PARAMETERS 

C " RETRIEVE COST FUNCTION DATA FROM ACSL RUN 
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i TWO WAYS TO CHANGE ACSI, .ARAHETERS PRIOR TO A 

r .^i^r^ nMTY FOR 1ST RUN AND SENSITIVITY RUNS 

C PARAMETER NEWRUN 



Ot NO Z2C0H ""tLE SAVED 

22COM SAVED INTO TMPFILB 



^S^lSir'^A^Ws-FOR- 
INTEGER NEWTON 

CHARACTER*10 FILE mrpmaX/1600/ 
INTEGER ICORRECT/1/, HBPMAX/lbuu/ 
LOGlCAl, CLHET/ .FA1.se. / 
REAI. TSTP/-35./ 

cIIIII- ZZCOM'FlirhSnaIIng"prior_to_^ 

c~IIirRUNF^Glir'IIlo3'us« parameter change, directly in ACSL 

IF (RUNFLG.EQ.l) THEN 

CLMET = -FALSE. 
END IF 

^ RUNFLG«2s Restore TMPFILE file 

IF fRUNFLG.EQ.2) THEN 

C.™-^ RUNFLG.3: Res^re ICFILE file, allow temporary ACSI, changes 



IF (RUNFLG.EQ.3) THEN 
IF (CLMET) THEN _ 
FILE « THPFILE 
TSTP » 200. 
CLMET » .FALSE- 

^^^FILE « ICFILE 

END IF PTTFi 8 Equ to ACSL> RESTOR 'ICFILE 

CALL 22SVRS2(0r FILE I 
END IF 
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RUNFLG-4: Restore TMPFILE file, allow temporary ACSL changes 



IF (RUHFLG.EQ.4) THEN 
CLMET » .TRUE. 

C IF f^f;;?"^^Q:^|J«^^nter TSTP for closed-loop runs:', TSTP 
C READfSr*) TSTP 

C END IF 

CAi:^ Z2SVRS2(0. FILB> « Equ to ACSL> RESTOR 'TMPFILE' 

END IF 

e::::::"£oarMEr;aram"in;rri5OTr^^ 

CALL SET_ACSL PARAM ( 



I 
1 r 
Ir 
Ir 
Ir 
1. 
I, 



RUNFLG 
VWXTAB 
VWXHO 
VWYAO 
KXER 
KAYER 
KRRER 



) 



RADPLG 
VW-iTAB 
VWXAO 
WVM I. 
KVXER 
K2ER 



RNSEED 
LVX 
VWXHl 
VWYAl 
KAXER 
KVZER 



NLAYER 
LVY 
VWXAl 
CLMET 
KYER 
KAZER 



RHOTAB 
LV2 
VWYHO 
TSTP 
KVYER 
KRGER 



IF ( RUNFLG. NE. 2) C^L 22EXEC 



I ACSL user interface 



IF ( RUNFLG. EQ.l) THEN 

CALL Z2SVRS2(1, ICFIL5> 
CALL 22SVRS2(le ICOFILE) 
CALL 22SVRS2(1, TMPFILE) 

END IF 

CALL 22SIML 



I Equ to ACSL> SAVE 'ICFILE' 

Equ to ACSL> SAVE 'ICOFILE; 
! Equ to ACSL> SAVE 'TMPFILE' 

I Run ACSL dynamic simulation 



CALL GET.AC«^PARM1( ^^^^ , COSTVY , COSTVZ , NCOST 

1. KXER , KVXER . KAXER . KYER , KVTER 

1 KAYBR r KZER , KVZER , KftZER , KRGER 

I KRRER RHOTAB , VWXTAB , VWYTAB ) 



IF (ROHFLG.EQ.l .OR. RUNFLG. EQ. 3) CALL 2ZBXBC I ACSL user inter 



cIIIIII ZZCOM File handling after completion^of^ACSL^run^ 



IF (ROHFLG.EQ.l) THEN 



CALL SAVENOMINAL( ICFILE » ^ 
! Make TMPFILE identical to ICFILE 



Save all xxxNOM arrays into ZZCOM 
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CALL Z2SVKS2(0r IC^^tf 
CALL SZSVRS2(1. TMPFILC 
BHD IF 



KVXER 
PHOTAB 



CALL '^^T.XC^PAR^M ( 

. KXER 
1* KAYER 
' KRRER « 

ELSE'ir fMEWRUN.EQ.O) THEN 
RETURN 

^^^CALL ZZSVRS2(lr TMPFXLF) 
END IF 

RETURN 
END 

SUBROUTINE GET_USER.PARAM 



NBPMAX . NEWRUN ) 



COSTVY 
KAXER 
KV2ER 
VHXTAB 



COSTVZ 
KYER 
KAZER 
VWTTAB 



NCOST 
KVYER 
KRGER 



C REQUESTS THE FOLLOWING P' 
OPTION 

RADFLGr RNSEED 
RHOUL. RHOLL 
VWXULr VWXLL 
\WYUL, VWVLL 
VWXHO, VWXAO 
VWXHi, VWXAl 
VWYHO, VWYAO 
VWYHl, VWYAl 



IN PARAMETERS: 



™« li Wi Ml 

VWY LIMITS FOR.|^S5_i:S?^** 
WIND PROFILE PARAMETERS 
WIND PROFILE PARAMETERS 
WIND PROFILE PARAMETERS 
SlND PROFILE PARAMETERS 



n 

Q 

c 
c 



lAYER ^ 
LVXr LVY, LVZ 



^11 hIIIhts in cost fcts 



IMPLICIT NONE 

INCLUDE •TAR.COHIB .FOR 

S;}' SliiSS: giS?S2; SSIilS 

IF ( OPTION. EQ.O) STOP 
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C IF f RUNFLG.EQ. 1) THEN roSTVZTHSL 
C PRINT 5, LAYER, LVX, LV2, C0STVZTH5L 



> ReA0<5,')' layer: Vii, '>y; LV2, COSTVZTHSL 
C END IF 

I rORMAT( /IX, 'OPTION . ' , 12 , 5X , ' RADFLG « ' . 12 , 3X , ' RNSEED 

5 ^' ^'F6RkAT( /IX, 'LAYER = ' , F5 . 0, 2X , 'LVX , LVY, LVZ , C0STV2THSL 
1, 3(F5.1,2X),F10.5, 'T') 

C RHO Limits 



IF 



(RUMFLG.EQ.1> THEN 
DO I 1. 40 
RHOTAB(I) 

END DO 
END IF 



RHOTABO ( I ) 



PRINT *r 
PRINT 10, 
PRINT 10, 
PRINT 10, 
PRINT 10, 
READ ( 5 , ♦ ) 
READ ( 5 , • ) 

10 



Enter RHOUL and RHOLL . 
'RHOUL (RHOUL (I),I«1,20 

'RHOTABO 25S5^0}f|'?:liS 
'RHOTAB « . (RHOTAB (I)»I*i»52} 
'RHOLL (RHOLL U).I-ir20i 

(RHOUL (I) 'I'f '20{ 
(RHOLL U)'^»^'20) 



I Upper limits 

I True RHO 

I Current Estimate 
! Lo%fe.r limits 



FORMAT{lX,All,3X,l0(F5.2,XX),:,2(/15X,10(r5-2,lX))) 



Q . VWXr VWY Nominal anri Limits 



PRINT 15 
READ( 5 , * ) 
PRINT 16 
READ{ 5, • ) 
CALL GET_ 



VWXHOr 

VWXHO, 
VWYHO, 
VWYHO 



VWXAO, 
VWXAO, 
VWYAO, 
VWYAO 



VWXHl, 
VWXHl, 
VWYHl, 

VWYHlr 



VWXAl 
VWXAl 
VWYAl 
VWYAl 



WIND TABLE ( RUNFLG- 1 ) 



15 
16 



FORMATC /IX, 
FORHAT( /IX, 



(^WXHO 
I VWYHO 



VWXAO) , 
VWYAO ) 



(VWXHl 
(VWYHl 



VWXAl) 
VWYAl) 



2(F10 
2|F10 



.0,F7.2) 
.0,F7.2) 



PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
READ(5, 



C 

C PRINT 
C PRINT 



10, 
10, 
10, 
10, 
) 



C READ(5,*) 



Enter VWXUL and VWXLL ?' ^ , 
'VWXUL (VWXUL (1)^1=^^20 

'VWXTABO VWXTAB0(I),I-lr20 
'VWXTAB ='r (VWXTAB (I)fl*l'20 
'VWXLL (VWXLL Ui,I-lr20) 

(VWXUL (D.I'^l'^O) 
(VWXLL (I), 1=1, 20) 



PRINT Enter VWYUL and VWYLL ?' 

10 'VWYUL (VWYUL (D'l-i'iS 

io' 'VWYTABO VWYTABO I)rl-lr20) 



! Upper limits 
! True vwx 
I Current Estimate 
t Lower limits 



i Upper limits 
! True VWY 
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PRINT 10 f 
PRINT 10 r 
PEAD(5r*) 
READ(5r*) 

PETURN 
END 



•VWYTAB 

'VWYLL 

(VWYUL 



(I), 1=1, 20) 
(I).t=l,20) 



(VWYLL (I)rl«lr20) 



i Current Estimate 
I Lower limits 



SUBROUTINE GET.WINO^TABL ? { K ) 



SUBROUTIHE TO COMPUTE BREAKPOINTS FOR VWXTAB AHD VWYTAB, 



C 

c 

r 

Q 

n 
C 

Q 

c 
c 

r 
C 
C 



BASED OH 



- lAYER 

- VWXHO, VWXAO 

- VWXHl, VWYHl 

USING polynomial CURVE StiAPB. 

XF K.Oe WAD VWXTAB INTO VWXTABO AND VWYTAB INTO VWYTABO, 



iKF'T^!cOMl8.FOR' 
REAL ALT, DAr DH 



C- 



— VWXTAB 
DA « VWXAl - VWXAO 
DO I*lr20 



20 ' ^°*='P ®" Breakpoints 

* I Altitude at breakpoint 

ALT « lAYERMl-l^^^^j^^^j^/^^xHl.VWXHO) 

IF (ALT.LT.VWXHO) THEN , constant wind below HO 

VWXTAB( X > ■VWXAO 

ELSE IF (ALT. GT.VWXHl) THEN , Constant wind above HI 

^ VWXTAB { I ) »VWXA1 



ELSE 



'^XTABdl-VWXAO - 3*DAMOH)"2 • 2*DAMDH)"3 
END IF 



IF (K.EQ.O) VHXTABO(I) « VHXTAB(I) ^ True wind parameters 



END DO 



VWYTAB 
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DA = VWYAl - VWYAO 

00 1=1,20 ^ ^^^P Breakpoints 

I Altitude at breakpoint 

IF (ALT.LT.vwYHO) THE" , Constant wind below HO 

VWYTAB ( I ) »VW> AO 

ELSE IF (ALT.GT.ywYHl> THEN , Constant wind above HI 

VWYTAB ( I ) • V. A 1 
^"^YTAB(I).VWYAO - 3-DAMDH)"2 • 2*DAMDH)**3 
END IF . 

IF (K-EQ.O) VWVTABOd) = VWYTAB(I) t True wxnd 

END DO 

RETURN 
END 

SUBROUTINE PR!^T_COST(K> • 
I SUBROUTINE TO PRINT REST.ir.TS AT LAYER ^KLAYER^ 



IMPLICIT NON*^ 

INCLUDE 'TAR_COMl8.FOR' 



I « KLAYER 



I Index to table 



RHOPCT 
VWXPCT 
VWYPCT 



, T . - / PWOTAB r I \ -RHOTABO ( I ) ) / RHOTABO ( I ) 
[l\ ' WXTAB I -VWXTABO I) 
\\\ Z l^;SYTABlit-VWYTAB0}li 

GOTO (100, 200, 30 J r 400) K 

ir TT AVER RHOPCT { I) r SQRT ( COSTVZ ( ILAYER ) ) 

RETURN 

««« TTKVvn VWXPCT(I)r SQRT (COSTVX( ILAYER ) ) 

RETURN 

300 PRINT 991 > IW^YER, VWYPCT(I), SQRT (COSTVY( ILAYER) ) 
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WPITB(40.991). K, IIAYEP. VWtPCT(I), SQRT( COSTVY( IIAYER) ) 
RETURN 

400 PRINT 995 . SAoPCTm , SQRTCCOS^ 

}' VWXPCT(I), SQRT(COSTVX IIAYER ) 

1 e y„y p,jx j I ) , SORT ( COSTVY ( llAYER ) ) 

1 t 



091 ^«^^'^T(I2,:4,4X.2rj0.3.;;;) 
QQ5 FORMAT(X2rl4r4X,6F10.3r r } 



995 

RETURN 
END 

SUBROUTINE SET METS 

§*iHirst.EFOUTIHE AliOWS C«ANGB TO MET PARAMETERS FOR COMPARISON 
C RUN TO NOMINAL 

**************************************************** 

TMPLICIT NONE 

INCLUDE 'TAR COH18.FOR' 

INTEGER lANS" 

PRINT • MODIFY WIND PARAM (0) OR TAB (1) 

READ(5r*) lANS 

C Wind 

IF (lANS.EQ.l) GOTO 100 
1 FORMAT(/A45,2(F8.1,2X,r5.2)) 

PRINT I, ' TAR> MODIFV ( VWXHO VWXAO ) , (VWXHl VWXAl) ? ' 
^ iT VWXHO. VWXAO, VWXHlr VWXAl 
READ(5,*) VWXHOr VWXAO, VWXHl, VWXAl 

PRINT 1, ' TAR> MODIFY fVWYHO VWYAO), (VWYHl VWYAl) ? 
^^l, VWYHO, vStAO. VWYHi. VWYAl 

READ(5,*) VWTHO. VWYAO. VWYHl, VWYAl 

CALL GET_WIND_TABLB(l) 
GOTO 200 

C — 
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2 rORMAT(/A20,l0F5.2,:.2(/20X,10(r5.2,lX))) 
,00 PRINT 2, 'TAF> VWX -T (VWXTAB, I) , 1-1,20, 

READ(5,*) (VWXTAB(I),I-1.20) 

rPINT 2 • TAR> VWY = ? ' , (VWrTABC 1 ) , I-l . 20) 
(- Change RHO 

'200 PRIHT 2. . TAP-- «H0 - ?',(WiOTAB(l,.X-l,20, 

PEAD(5,M (RHOTAB(l). 1=^.20) 

nFTURM 

^ REAL FUNCTION FLNEAR( NBPMAX, NBP , TAB , IHP ) 

r HBPMAX s waxiwuni immtaer of breakpoints 

r NBP '• current number o£ breakpoints 

C TAB : table oLSSy?*^^®^" - values of dependent variable 

I ?^SiipMA^'i2-NBPMAX) ; vStSIs of inSependent variable 

\ \ii outU"!rtH2^3ipiiidtsrra?ti5i^^^^^^ " --^ 

IMPLICIT NONE 
REAL TABCK^INP 
INTEGER NBPMAX, NBP, I 



C interpolate/extrapolate using the first slope 



FLNEAR - TAB^xj ^TAB(NBPMAX*2)-TAB( 



- (TAB<2)-TAB(1) ) 
NBPMAX^ 1) ) 



RETURN 
END IF 



C Find index of independent variable value greater than .nput value 

°° if' {INP.lI!tAB( NBPMAX*!)) GOTO 100 
END DO 

C interpolate (extrapolation ueing last slope for input values 
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C larger than last breakpoint value) 

100 CONTINUE ^ ,TNP-TABfNBPMAX+l-l)) * (TAB( I ) -TAB ( I-l ) ) 

FLNBAR • TAB(I-X) *. ' JtAbThBPMAX*!) -TAB(NBPMWt+I-l ) ) 

RETURN 



^ SUBROUTINE LSQ_P0I.Y2 ( C - A , B , N ) 

C points (Air Bi) • 
C 

C Outputs 

r 

C Inputs s 

C A(N) : vector of values of Independent variable 
r B(N] : vector of values of dependent variable 
r N ! number of data points 

j*********.******^****-'^**** *********** 

IMPLICIT NONE 

REAL AMAT(3,3), BMAT(3). AMAT.INV( 3 , 3 ) , DET.AMAT 
INTEGER Ir K 
r— - — - — - — 

C — Compute AMAT — — 

r • compute 1st column of AHAT 

DO K =^ I r 3 

AMATfK,!) « 0.0 
DO I « 1 r W 
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AMAT(K,1) -AMAT(K,1) -A(I)**(5-K) 
Zm DO ! I 
END DO ! K 

C Compute 2nd column -^f AKAT 

AKAT(1,2) - AMAT(2a) 
AMAT{2,2) « AMAT(3,1) 
AMAT(3r2) « 0.0 

AMAT(3:27 - AMAT(3.2) - A(I) 
END DO ! 1 
C — - Compute last column of AMAT 

AMAT(1.3) = AMAT(2,2) 
AMAT(2,3) = AMAT(3,2) 
AMAT(3,3) = N 

• " — 

C Compute BMAT — 

C — * 

BHAT{ 1) = 0.0 
BMAT(2) =0.0 
BMATC 3) « 0.0 

BMAtIi) ^ BMAT(l) * B(I)*A(I)**2 
BMAT 2 » BMAT 2 ^ B(T)*A(I) 
BMAT(3) • BMAT(3) * B(I) 
END DO ! I 

r — 

C Invert AMAT ^ 

CALL INVERT3BY3( AMAT , AMAT^INV , DET^AMAT ) 
cI"ZII*Compute polynomial coefficients^ 

DO I = 1 , 3 
C(I) « 0.0 

C(I)^»'C(I) * AMAT^TNV(I,K)*BMAT(K) 
END DO ! K 
END DO i I 

RETURN 
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c-BROUTINE INVERT3BY3( A , A_INV , DET.A ) 

? Tliis routine inverts a ^^iJ^^^T " , ^ 

I z TaH: l^l 5l^^lnSnr.ll SeT.A is set to zero if A is almost 

" singular. 



C 



IMPLICIT NONE 

REAL M3,3), A I''X^nET*^ 

^u^ii'JT: iFifSoS/ 1.0.20/ 

INTEGER Ir J 
Transfer matrix A into work array 



DO I « 1 r 3 

END no 
END DO 



DET_A « D_DET_A 
Q Test for singularity 

IF (ABS(D_OET A) .XE .EPSILON ) THEN 

DET A = O.U 

RETURN 
END IF 

C compute the matrix of cofactors of A 

i§?{3' 1 : ^ g-AU:l!*D:A|^Ji - DIa12,2{*D.aIu3) 
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Compute the matrix Inverse of A 

00 I » 1 * 3 

A INV(iJ) » COr(J,l)'0_DET.A 
END DO 
END DO 



PETURN 



C 
C 

c 



SUBROUTINE PL0TSM0OTHED( SMOOTHED, X, N) 



Inputs! 



C SMOOTHED 
C X 
C Y 



smoothed dependent variable (dim N) 
jSSSSSSrvariSfe?^'§eiSi2 Smoothing (dim N) 



INCLUDE 'COMPLOT* 

INTEGER N 

REAL SMOOTHED! ♦Ir 

DATA DATLEN/lfiOO/ 



NAME(115) - 'XSMOOTH' 
NAME 116 - 'YSMOOTH' 
NAME(117) « 'YNONSM- 

DO I « 1 r N 

CALL SAVEDATA( 
CALL SAVEDATA( 
CALL SAVEDATA( Y(I) 

END DO 1 I 

CALL PLOTROUTINE 



X(M. Y(*) 



X(I) 

SMOOTHED ( I ) 



DATLEN , 
DATLEN r 
DATLEN, 



I. 



115) 
116) 
117 
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PFTURN 
END 



- ^TiBPOtrriNE HOV^AVEHAGEf Vi, YO, N, WP, P) 
0 This routine compute i:He sequence Yl as follows: 

'r vnci), i-UN : input sentience 
r wf 1 ) . j«i.M s wexgnt? 

C Yl(l) 
C Vl(N) 

C Yl(i) « ( Y0( 
r [ WP(] 

? - • 



IMPLICIT NONE 

INTEGER I* J 
REAL WO, Wl 

C Initialize output sequence 

DO I = 1 r N ^ 
Y1(I) - 0.0 
END DO 

Q Compute sum of weights 

wo « 0.0 

wp'jf - ABS(WPfJ)) 

= wo + wp(j) 

IF^(WO.EQ.O.O) RETURN 

r — Compute output sequence Yl 

Yl(l) « 0.0 
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Y1(N) = 0.0 
Wl » 0.0 

YlTll » Yl(l) * Y0(J>*WP(J) 
Yl N = Yl N * YO(N..)-i)*WP(J) 
v:i « Wl WP(J) 
END DO 

IF (Wl.NE.0-0) THEN 

Yl{i) « Yl(l) / Ml 

YUN) « Y1(N) / Wl 
EHD IF 

DO I * 2 , N-1 

Y1(I) = Y0(I)*WP(1) 

Wl • Wl * WP(J> 
END DO 

DO J - 2 , MIN(P,N-I-^l) 

iMD « Y1(I) + YOi-t-^-J-D^WPt J) 

Wl - Wl * WP(J) 

ir^(Wl.NE.O.O) Y1{1) = Y1{I) / Wl 
END DO 

RETURN 



c 

f ' 



SUBROUTINE RADARSET( PROJ^DIAM, KSEEDr SETFIAG) 



C This routine calculates: 
C - the range resolution, , ^. ^ pp*r 

r - «"ho T-jidlai veiocitv resolution due to rtT, , 

C - and "e free space S/fi at a standartd range of one nautical mxle. 

INCLUDE 'RADAR COM18.FOR' 
INTEGER KSEED.^SETFLAG 
REAL TRF, PRF . DF 
REAL PROJ.DIAM, AG 

Q Store argument values into common block 

PROJ DIAMETER - PROJ^DIAM 
ISEED » KSEED 
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C Ask user radar settings 

IF (SETFLAG-EQ.l) THEN 



PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT ♦ 
PRINT ♦ 
PRINT * 

READ(Sr 

• r 



WEATHER 0 (Clear) 1 

WAVE LENGTH (m) 

RECEIVER y^nZQ {Hz) 

RADAR POWER (W) 

ANTENNA GAIN (dB) 

NOISE FIGTjnC (dB 

TOTAI* LOSSES (dB) 
FFT LENGTH 



(rain) : 



WEATHER 
WAVE LENGTH 
RECETVER^FREQ 
RADAR POWER 

antenfa gain 

NOISE ETGURE 
T0TAL"L0SSES 

fft length 



TOTAI^lLOSSESr FFT.LENGTH 



SETFLAG « 0 
END IF 

C compu..e beamwidths assuming they are equal in azimuth and elevation 



AG » lO.**(0.1*AirrBHNA GAIN) 
ELBV^BEAHWIDTH " A2IM_BEAMWIDTH 



OEG2RAD 



Compute range resolution 
RANGE RES - 150. B6 / RECEIVER.FREQ 
C "compute radial velocity (range rate, resolution <£or 20% duty cycle 

TRF = 5.E6 / RECSIVER.FREQ 

PRF « i - B6 / TRF 

DF = PRF / FFT_U5NGTH 

RRATE^RES - 0.5 * WAVE_I.ENGTH • DF 

C PRINT *, ' DF OF, • DV ' , RRATB.RES 

C Compute maximum radar cro»» section 

liiSi^" : §:f *%!^%Ic?f.5^^rgX;^_«NGTH"2 

RCSMAXDB = 10. ♦ ALOGIO ( RCSMAX ) 
C compute free space S/N ratio at 1 nautical mile range for maximum R 



STN^O - 10. : A^O<=iO(RAOAR^POH|R,^^^ 
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- 10. * ALOC-U»(RECEIVER FREQ*l.E-6) - 60. - NOISE.FIGURE 
: TOTAL_£oSSES * 10. • XLOGIO ( FrT.tENGTH ) 

r PRINT * , • STN_0 - STN_0 1 ! 

r Save values for display 

^"^'^ ^WeItHEr! " WAVE r.ENGTH, RECEIVER FREQ 

" oSnlS POWER ANTEWNA GAIN. NOISE FIGURE 

" TOTA?-£osiE3. FFT LENGTH, ELEV^BEAHWIDTH 
; [ AZIM_5EAMWIDTH ) 

CALL INFONAMES 
r Identify frequency band 

ELSE IF t„^OrEBNGTS.GT.LIGHTSPEE^ 

PRINT BAND, wavelength » ' ,^Ave_lewgth 

^^^PRINT ♦ 'FREQUENCY BAND NOT IN DATA BASE FOR' 

PRINT 'WAVELENGTH MAVE.LENGTH 
END IF 

Q Compute* attenuation due to the weather 

i^^t™V^El^.^lH^"^LIGHTSPEEO/lCA_BANO(l)) THEN ^ ^^^^ 
EI,SE"PywIvE\ENlTS.LE.LIGHTSPEEO/K_BAMD( 1) ) ,THEN^^ 
ELSeIP ( hIvE^'lENCTH . . LIGHTSPEED/KU.BAHD ( I ) j THEN 
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ATTENUATION « 0.4 

ELSE « , 

ATTENUATION » 0-1 

END IF 

^^^ATTENOATION - 0.0 
END TF 



PRINT * 
RETURN 



'ATTENUATION * 



! KU BAND 

I X BAND AND BELOW 



ATTENUATION 



" SOBROUTIHE «ADAR«EA|( |^DAR FLAG, ASPECT^ANOLE , HM 

:: rInoiIe. rSioiIb, BUIOISE, AZMOISB) 

c ' •••••••••••♦•♦♦♦•♦•*•*•*••••••**' 

INCLUDE 'RADAR COMIS-FOR" 

11^ ^?'hI™I;NIKV. STN, GAUSS. RCS.DB, STN^DB, KSIGMA 

'"'^'^ 2^§F'^iANGi''lA?Ef E^^^ AZIMUTH 
SJHJpfi RRMEA? ELMEAS, A2MEAS 

S^SoiIe, IrNOISE. X2NOISE, ELNOISE 

C Test for RADAR^FLAG 

IF (RADAR.FIAG.BQ.O.O) THEN 

HM = 1- 

RGNOISE * 0.0 

RRNOISE « 0.0 

ELNOISE « 0.0 

A2NOISE « 2;S,-«. 

RGMEAS « RANGE 

RRMEAS - FANGE^RATE 

BLMEAS - ELEVATION 

A2MEAS - AZIMUTH 

RETURN 
END IF 
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r Compute the normalized radar cross section, RCS_NORM 

ANG - PI • PROJ DIAMBTCP / WAVE LENGTH 
ANG ri r-nw-^ s IN ( ASPECT_ANCI:E ) * l.E-5 

P.;S NOBM = ABS( SlN(ANg)/\NG ) * 0 01 
RCS'DB - 10. •ALOGlO<RC?_l«JRM) 

C PRINT *. • RCS_NORM RCS_NORM It 

L- Compute attenuatici' clue to the weather 

J: &"ved to routine RADARSET and Kioete 

r .... viriable ATTENUATION is passed via a cononon block 



I l^^tSAVrEgNliH^^f'l.. THEN 

I ELsf IP^ywI^^ENlTS . UE . 1 . 5 ) THEN 

I ELsl^rVwI^^ENlTi; . LE . 3 . , THEN 
C ATTENUATIOfl « 0.4 

r ELSE 

C ATTENUATION « 0 . 1 

C END IF 

C ELSE ^ ^ 

r ATTENUATION -0.0 

C END IP 

r Compute the actual S/N ratio 

STNDB = STNF - l.E-J * ATTENUATION * RANGE 
STN_DB» STNDB 

C PRINT *r • STNDB * r STNDB 

C Compute the p:robability of detection 

IF {STNDB.LE.4 0 THEN 

PROS DETECTION « 
ELSE IF"(STNDB.LE.6. > THEN 

PROB^DETECTION « 6.1 
ELSE IF (STNDB. LB. 8. > THEN 

PROB DETECTION • 3 
ELSE IF"(STNDB.LE.10. )^THEN 

PROB DETECTION « 0-6 5 
ELSE IF""( STNDB. LE. 12. ) THEN 
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ELSE IF"~(STNDB.LE.14.) THEN 
PROB.DBTECTION - 0.92 

^^^PBOB DETECTION - 0-?9 
END IF ~ 

set randomly the H«.t/Miss flag 

RH - RAH(ISEED) 

IF fPROB DETBCTIOH. GT.RN) THBH 
HM - T.O 

HH » 0.001 
END IF 

^ If RADAR_FIAG-1. cn.ly HM flag is used (no noise) 

IF (RADAR FLAG.EQ.l.) THEN 

RGMEAS"- RANGE 

RRMEAS - RANGE RATE 

EI.MEAS - ELEVATION 

AZMEAS - AZIMUTH 

RETURN 
END IF 

e/M to 30 dB (to limit the system accuracy 

E:::::: ^|"i%^S!iSnr?!a!;le,'!S^%h^ rUe ^ate signal 

e compute the error .n range rate measurement due to g.ussian noise 

^ I uniform dxst. 

tJM « RANf ISEED) - 

luir-4QiT?:iN^TR^i>"'* C0S(THOPI-RH, « ^—i- dist. 

I llo^f l^i^G^r HfS(Kiiill^!Al^(^Sis) , . GAUSS , 
RRNOISE - GAUSS • RBATE.RBS • SQNINV 

5--:::: J2„SS.'S.ii?&T„,!rr. Sr.'' 

C angle signals 
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IF (STHDB.GT.40.) STNDB = 40. 

STN » 10.**(0.1*STHDB) 
SQNINV - 0.5 / SQRT(STN) 

compute the error on range «e.»ur«ni«nt due to g.ueaian noise 

i uniform dist. 



^lu^l^'-ioRif^iNSSGlRHiT-i COS(TWOPI*RH) « g««"i«" '^i"' 

g ^Ku^f !'i?S^r«Si^[iiiS^!^^S<GASsS,) , GAUSS , 
RGNOISE - GAUSS * RANGE_PCS * SQHIHV 

C compute the error azimuth measurement due to g.u.sian noise 

I unifoxin dist. 

RN = RAN(ISEED) ^ 

IKu^f =''i6Rif:2!*LSGlR^.M'* COS(THOPl-RH) « ^auasian diet. 

I lIu^fM^G^f^MfSil^lllSJ^lA^^^aASsS), , GAUSS , 
A2N0ISE - GAUSS • A2IM_BEAMWI0TH • SQHIHV 

c compute the error elevation measurement due to gaussian noise 

t uniform dist. 

^Kuls'^i^'ioRi'^iNSoGlRNiT'* COS(TWOPI-RH) • gaussian dist. 

g ^L^fMtSgrM?SiKiiiS^!ASt(GASsS)) . GAUSS ) 
ELNOISB - GAUSS • ELEV.BBAMWIDTH * SQHIHV 



set radar measurements 

RRMEA'S » RANGE * RGNOISE 

RRmIaS - RANGE RATE * RRNOISE 

ElS^ - ELEVATION + f^»;pl|| 

AZMEAS - AZIMUTH + AZNOISE 



RETURN 
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END ^^^^^^♦♦♦♦♦♦♦♦♦♦^ 

C* ************* 

*' SUBROUTINE IHFONAMES 

C ^^^^^w*^*****-****^**********' 
INCLUDE -COHPLOT' 

TNFOOl) • 'V»BATHER' 

IHFO 32 - 'WAVELEHGTH' 

INFH 3V) - 'RCVR FREQ- 

iwF" 34 - 'RAOARPOWER 

THF6(35 - 'AHTBN^GAIN- 

' 'fU'EIngth- 

INFO 3° = 'EL BMIOTH- 
INFO (40) = 'A2IBWIDTH' 

RETURN 

END -.*...*•♦••*•••••*♦••♦•♦•••• 



^' SUBROUTINE LSQ_POi;.Y( X. Y, W, B, PI, C lER) 
C least-square sense. 



c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

n 



Inputs 



X 
Y 
W 



4) N 

5) PI 

Output 



.... . B I iiSra !r,i 

SSSSr si ^SelflctSnli (i.e., polynomial order * 1). 



vector of 
vector of 



C 

lER 



vector of si^e PT containing the polynomial coefficients, 
IMSL error parameter 

129 - S2t?ix°A algorithmically singular 
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34 - accuracy test failed (number ofaccurate digits IDGT 
C not achieved) 

I internal variable, needed for the_ll«L_routine_LEOTlF 

8 rhrmir;;u;i;rsoive;-the e^«^^ ^« 

^ ii 11x8;"" ' vlct« of rlghlinand side values 

^ i I nSSber"? accurate digits for the value, of A and B. 

^ (hardcoded as 3) 



IMPLICIT NOHE 



• ^^lE^lW: ?tH). W(N), C(PI) 

REAL A(in,10)r B(10) 
INTEGEP Ir ' 
REAL WKARBA( 10 ) r XX 
INTEGER IDGT, lER 

C compute the elements of nhe first: column of A 



DO K 1 r PI 
A(K,1) « 0.0 
DO I - 1 r H 

iF^JtlQTo .AND. X(I).EQ.O.) THEN 

XX « 1.0 
ELSE 

XX = X(I)**J 
END IF , ,^,Tv*vv 

A(K.i) » A{K.l) * W(I)*XX 
END DO 
END DO 



C- 



compute the elements of the last column of A 



DO K - 1 f ^ 
A(K.Pl) - 0.0 
DO I * 1 f N 

If"("eQ.O .AND. X(I).EQ.O.) THEN 

XX - 1.0 
ELSE 

XX - X(I)**J 

A(kJx) - A<K,P1) - W(I)-XX 
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END DO 
END DO 

compute the upper left elements of A 

n'» t- = 2 r ^ 

A.(r.t.l « A(K1#1) 
END DO 
END DO 

r compute the lower vight elements of A 

no K * ri*2-L r PI 
■ r2 « K+L-Pl 

AfKcLI - A(K2,P1) 
END DO 
END DO 

rompute the vector H 



cc 

UO K = 1 r PI 

UO I = ' N 

'lF'("iQ.O .AND. K(I).EQ.O.) THEN 
XX » 1-0 

« x(i)-j 

|(K,'^ B(K) - Y(I)*W(I)-XX 

END DO 
END DO 

Solve for C: A*C « B 

IDGT * 3 

'I Bi in B. IDGT, WKAREA, lER) 
CAI^I- LEQT1F( Ar Ir Plr 10' 

C(I) » B(I) 
END UO 

RETURN 
END 
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TAR 8IMX8.CSL . 

^.Ttrr^BV CTMUUVTTON TAR SIM18.CSL 02-MAY-X990 

PROGRAM TRAJECTORY SIMUI*AriUN - 

Fortran progreun. 
Modifications s 

:- ^'■""P^o" Requires the macro file NONE 

AJLJ. velocities are in metres per second ^ 

or rad/ sec 

p»*oi-nal Units*. All angles are in degrees 

External units. «. listances in metres ^ . 

All velocities are in metres per second _ 
or deg/sec 

Source^ :^ Spacecraf t^viO , ^n6^June^l973^^pp384;ee • 

''°"-^;HJ4purjiiiH4r on input line 28 of page 27 of file -ANALYSIS: rARMAl..T 



M;;E;5"FKBirKB;:";i^raT"K5^C. KDAI, K0A2, KDA3, KDA4, KDA5) 

■ ?S!I^S;S;ro'com?«irthe aerodynamic coefficient KDA as a polynomial 
- function of Mach number. 

" Inputs: 

" : iJS^C ; X?«y"ST''"eakpoints «or Mach number 

• - KDA1-K0A5 t Arrays c€ polynomial coefficients 

• -^KDA*' ' Aerodynamic coefficient 



SUBSTITUTE SHEET 



PCr/US92/03377 

WO 92/19928 



62 



) ) 



iu^CRO RELABEL LI. " . L3 , t-» . LEND 
PROCEDURAL ( KDA ■ MACH ) 

IF I MACH .GT. KDAMAC(l) ) GOTO LI 

KDA - KDAI(l) 

GOTO LEND 

- U^-^i.rt-iiVi Ell! : iSS : toiStii'i i 

GOTO LEND 

||-(-^'^H"ViT KDA^tAC,^ GOTO L3 ^ . ^ ^oA3(3) - 

KDA - KDA3(1) I ^ j + MACH * KDA3(S) ) 

GOTO LEND 

^^f^^H^^GT. KDAMAC(4)J gTO L4 ^ . , ^^^^^^^ ... 

KDA = KDA4(1) t I^^^H ^ ( ^^4(2, ^ ^ KOA4(5) ) ) ) 

GOTO LEND 

"'K^Ji'^'iSisd) * KDA5{2)*«ACH 
LEND. .COKTIHUE 

END $ " of PROCEDURAL " 

MACRO END 

^'^°"-^?I^SrSSntro7inrut: line 78 of page 27 of file "ANALYSIS: (ARMAL.T 

■ ^Sil'SS4ro'cSSpG«Pt:he aerodynamic coefficient KDO as a polynomial : 
• function of Mach number. - 

- ^"SaCH* "Sf^'breakDoints for Mach number I 

; : lllft^Ol \ mih ot^lllylo^^^' coefficients . 

" S^kEo*" : Aerodynamic coefficient I 

MACRO RELABEL LI. I|2 . L3 , L4. L5, L6, LEND 
PROCEDURAL ( KDO » MACH ) 

IF f MACH .GT. KDOMACtH ) GOTO LI 
KDO » KDOl(l) 
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GOTO LEND 

if^i^ffi*""!!. KD0K|^(2)J GOTO L2 ^ . , ^^02(3) .•• 
KDO = KD02(1) : U^^U : S ^lllW) * MACH - KD02(5) ) ) ) 
GOTO LEND 

^l i^ffi^'^iT. KD0HAC<3,J GOTO L3 ^ ^^03(3) ... 

KDO - K003(i) : IJI^CH ^ { Kggljll ^ * kD03(5) ) ) ) 

GOTO LEND 

^^(•^?S''"iT. KrOHAC,4,J GOTO L4 ^ ^„ . , ^^04 ( 3 ) ... 
KDO - KD04(1) * HACH ^ | Kgg* | ^ ) ^ , ^004(5) ) ) ) 

GOTO I.EHD 
GOTO LEW 
GOTO LEND 

L6.. CONTINUE ^ KD07r2) MACH * ( KD07(3) ... 

KDO . KD07(1) * KACH ^ j ^ * KD07(5) ) ) ) 

LEND. .CONTINUE 

END $ " of PROCEDURAL " 

MACRO END ^ 

''''":S^^I^IpurSsi'1^rorinp«t ime 135 of p*,* 27 fUe -analysis: CARMM. 

: ^^H^^2;ri';oS?«*«'t»>« aerodynamic coefficient KtO .s . polynomial ; 

- function of Mach number. - 

-"PsEac • AlrVSi^breakpeints for Mach number 

" : KLOltKLOS ; A^«y8 of polynomial coeffxcxents 

- 2"kEo^' « Aerodynamic coefficient 
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MACRO REIABEI. LI, tj2, L3, L4 , LEND 
PROCBDURAI.( KLO " MACH ) 

IF ( MACH .CT. KLOMACM) ) GOTO LI 

GOTO LEND 



l^f'^ffi;;;!!. KL0tt^(2)J GOTO L2 ^ . , ^^,02(3, . 

KLO - 10.02(1) * MACH ^ ( ^ • KL02(5) ) ) ) 

GOTO LEND 

^|;;^SACr"iT KLO^^^ GOTO L3 , . ^ ^,^03(3) ... 



^^i'^M^H^YIr. KLOHA^^ GOTO L4 ^ ^^H - ( KI,04{3) ... 

GOTO LEND 
b4 . - CONTINUE n c / ^ ♦MarH 

LEND. .CONTINUE 
END $ " of PROCEDURAL " 

MACRO END * 

••''^^^-^;'i^tpTSS?i^rin"Lput line 181 of pag. 27 of file "AKALTSIS. {ARMAL. 

~ ?^Il'SS4«'c5S?G«rt:he aerodynamic coefficient KM as a polynomial : 
- function of Mach number. 



' -"!MSc i S^Sy"Sf breakpoints for Mach number 

: KMl^KMS t ^ays of polynomial coefficients 

- Outguts ^ Aerodynamic coefficient 

^CRO RELABEL LI, L2, L3. LA, L5, L6 , L7 , L8, LEND 
PR0CED0RA1.( KM - MACH ) 

IF ( MACH .GT. KMMAC(l) ) GOTO LI 

KM - KMl(l) 

GOTO LEND 
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LI 
IF 



.CONTINUE 
( MACH .GT. 
KM - KM2(1) 



KMMACf2) 
+ HACH * 
* MACH * 



GOTO LEND 
L: . .CONTINUE 
IF ( MACH .GT. 

KM » KM3(X) 

GOTO LEND 
L3. .CONTINUE 
IF ( MACH .GT. 

KM = KM4<1) 

GOTO LEND 
L4 . . CONTINUE 
IF ( MACH .GT. 

KM = KM5(i) 

GOTO LEND 
L5. .CONTINUE 
IF ( MACH .GT. 

KM - KM6(1) 

GOTO LEND 
L6 . .CONTINUE 
IF ( MACH .GT. 

KM - KM7(1) 

GOTO LEND 
L7 . . CONTINUE 
IF ( MACH .GT. 

KM « KM8(1) 

GOTO LEND 
L8. .CONTINUE 
KM « KM9(1) 

LEND. .CONTINUE 

END S " Of PROCEDURAl. 
MACRO END 



KMMAC ( 3 ) 

* MACH • 

* MACH ♦ 



KMMAC ( 4 ) 
+ MACH * 
MACH * 



KMMAC(5) 
* MACM • 
-»■ MACH ♦ 



KMMAC ( 6 ) 
* MACH ♦ 
+ MACH ♦ 



KMMAC ( 7 ) 
-t- MACH ♦ 
* MACH ♦ 



KMMAC ( 8 ) 
+ MACH * 
+ MACH * 



-I- MACH 
* MACH 



GOTO L2 
KM2(2) + MACH ♦ 
KM2(4) + MACH ♦ 

GOTO L3 
KM3(2) * MACH • 
KM3(4) MACH ♦ 



GOTO L4 
KM4(2) MACH * 
KM4(4) * MACH ♦ 



GOTO L5 
KM5(2) * MACH ♦ 
KM5}4) ♦ MACH • 



GOTO I#6 
KM6(2) * MACH ♦ 
KM6(4) + MACH * 



GOTO L7 
KM7<2) + MACH ♦ 
KM7(4 + MACH ♦ 



GOTO LB 
KM8(2) ^ MACH ♦ 
KM8(4) + MACH * 



KM9(2) + MACH * 
KM9(4) + MACH * 



( KM2(3) 
KM2(5) ) 



{ KM3J3) 
KM3<5) ) 



( KM4(3) 
KM4(5) ) 



( KM6(3) 
KM6(5) ) 



( KM7(3) 
KM7(5) ) 



( KMe(3) 
KM8 ( 5 ) ) 



) ) 



) ) 



) ) 



) ) 



) ) 



) ) 



) ) 

9(5) ) )') 



JNOFF-W.IIF, -L ignored 

on output page 38? on input 



line 248 of page 27 of file "ANALYSIS s (ARMAL. 



MACRO rKA( KA, MACH, KAMAC, KAl , KA2 ) 



^Ail'SS^ro^cSm^Gtes'the aerodynamic coefficient KA as a polynomial 
function of Mach number. 
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■ inputs: 

• - MACH 
. KAMAC 

• KAi-KA9 

• Ounput: 



A^2y"Srb.veakpoints for Mach number 
: A^eys of uolynomial coef fxcients 

t Aerodynamic coefficient 



MACRO REIABEL LI, LEND 
PROCEDURAL ( KA ■ MACH ) 



IF ( 
KA 



MACH .GT- KAMACrl) | GOTO LI ^ ^ 



GOTO LEND 

LI.. CONTINUE 



END 



KA = KA2(X) * JJ^LH * ( KA2|2) ^ n^^^ ^ kA2(5) } 
LEND. .CONTINUE 
- o€ PROCEDURAL - 



) ) 



) ) 



MACRO END — — 

-^^^^t^^^^o. incut line 281 of page 2. of file -AKALXSIS^ARMAL. 



^Sfir.J^ifSS • 

earth coordinates of wind velocity relative to the Earth ,n„sec, 
TABLE VWXTAB* 



0 

SQOO . 

i.0E4 

1.5E4 

7.07 

7.07 

7.07 

7.07 



1000. 
6000. 
1.1E4 
1.6E4 
7.07 
7.07 
7.07 
7.07 



1, 20 • 
, 2000. 
, 700f>. 
, 1.2E4 
; 1.7E4 
, 7.0" 
,7.07 

, 7.07 



TABLE VWYTAB* 1 



0. 
5000. 
1.CE4 
1.5E4 
7.0" 
7.0" 
T.0-? 
7.0T 



1000. 
6000. 
1.1E4 
1.6E4 
7.07 
7.07 
07 



7.07 



2 0 / 
2000- 
700'» . 
1. 2E4 
1 . 7P:d 
7 . 07 
7 .07 
7.07 
7.07 



3000. 

6000. 

1.3E4 

1.8E4 

7 .07 

7.07 

7 .07 

7.07 



3000. 
8000 . 
1.3E4 
l.eE4 
7 .07 
7 .07 
7 .07 
7.07 



4000. 
9000. 
1.4E4 
1.9E4 
7 .07 
7.07 
7.07 
7.07 



4000. 
9000. 
I.4E4 
.9E4 
.07 
.07 
.07 
.07 
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Air density (kg/in* 



2) 



TABLE RHOTAB, 1 



0. 

5000. . 
1.0E4 . 
1.5E4 . 
1.2245. 
0.7357. 
0.4X12, 



1000. 
6000, , 
1.1E4 , 
1.6E4 , 
1.1095, 
0.6559, 
0.3616* 



20 / 
200U. 
700n . 

1 . 71^4 
1.0031, 



3000. , 
8000. , 
1.3E4 , 
1.8E4 , 
0.9082, 



0.?f?72, 0.5218, 
0.3118, 0.2629, 



4000. 
9000. 
1 .4E4 
1.9E4 
0.8184 
0.4636 
0-2259 



o:i929: o:i64o: 0:^23: o.im; 0.1029/ 

Tempera'jure (degK) 



TABLE TMPTAB, 1^ 20 / 

0. , 1000. , 200^. , 

5000. , 6000. , 7000. , 

1.0E4 , 1.1E4 , 1.2E4 , 

288.20, 281.69, 215^20. 

255.70, 249.10, 242.63, 

223.25, 216.75, 215-05, 

216.70, 216.70, 216.70, 



3000. 
8000. 
1.3E4 
1.8E4 
268.70, 



4000. 
9000. 
1.4E4 
1.9E4 
262.20 



236.35, 229.70 
216.65, 216.70 
216.70, 216.70/ 



Correction for wind coordinates and air density 



TABLE VWXCOR, U 1600 
TABLE VWYCOR, 1, 1600 
TABLE RHOCOR, 1, 1600 



0., 
0., 
O.r 



1599*1. E6, IfOO^O. / 
1599*1. E6, 1600*0. / 
1599*1. E6, 1600*0. / 



'""""o^'outp^t pigrio; on input line 340 of page 27 of file "ANALYSIS: (ARMAL. 



CINTERVAL 

CONSTANT 
CONSTANT 
CONSTANT 



CINT =0.1 

DTORAD = 57.295779 
G = 9.80665 

PI « 3.1415926 



$ - BRL CONSTANT CI 



CONSTANT £fTfXc^«\. S -0.9627 BRL multiplier on KLO and KLA " 
COnItaSt CLPFAC - 1.0\ CXOFAC « 1.0 , CX2FAC « 1.0 

'^^"-S^'oStpi;^ j2Se"o; on input line 357 of page 27 of file "ANALYSIS: ( ARMAL. 
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TMITIAL 



BlUi Aero 



data from BRL FCI-155 * 



ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 
ARRAY 



KDAl ( 1 ) - 
KDOX(l) . 
KD07(5) , 
KL01(2) , 
KMl ( 1 ) r 
KM7 f 5 ) r 
KAl ( 5 } f 



KDA2(?), KOA3(5 
KD02(J), KD03(5) 
KO0MAC{6) 



KL02(5l 

KM2 ( 5 ) I 
KHO ( 5 ) 
KA2 ( 3 ) 



KL03(5) , 
KM3 ( 5 ) r 

KiCMAC( 1) 



KDA4 ( 5 ) r 
KD04(5) . 

KL04^5) , 



KDA5 ( 2 ) , 
KD05(5)< 

KI,05(2) 
KM5 ( 5 ) r 



KDAMAC(4 ) 
KD06 ( 5 ) 

KL0MAC<4) 
KM6(5) 



Define constiants in 



Aero Data Arrays 



CONSTANT 
CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 
CONSTANT 

CONSTANT 
CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 
CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 
CONSTANT 

CONSTANT 
CONSTANT 



KDAl = 

KDA2 = 

KDA3 = 

KDA4 « 

KDA5 = 
KDAMAC« 

KDOl « 

KD02 « 

KO03 * 

KD04 « 

KD05 ^ 

KD06 « 

KD07 = 

KDOMAC- 

KLOl - 

KL02 - 

KL03 = 

KL04 = 

KLOS « 
KLOHAC' 

KMl 
KM2 



-28 .981984 

100.48648 

1.89 . 37350 

-568-33625 

-37.317582 

31.613187 

5.6409091 

0.80 

1.80 

0.05 

106.55819 

-670.28411 

-621.41655 

3216.5548 

103-02763 

-109.02406 

-1580.9899 

5271.5295^ 

-5.2380076 

9 .4212207 

0.20694309 

0.034380058 

0.84 

1.035 

0.55 

5.2834534 
-62.964179 
157 .96430 
-747.59628 
-30.483617 
75. 152710 

: 0.71218182 

= 0.60 
1.23 

r 1.28 

t 1.7205997 



129.07274 

•18.128487 

-694.01933 

126.47310 

94.483221 

-4.6953080^ 

-0.90909091 

0.97 . 



-495.42480 r 
194.96965 
2716.0843 , 
-871.42812 
-315.74748 r 
0.0 

5951.9060 

-1240.2931 

16.492312 

-1.7685151 

-.0040499634, 

-.0046420421 

0.915 

1.08 ' 

0.057692308 
-28.181146 , 
23.626028 
-669.99668 , 
196.31506 
105.03552 
-15.991247 
0.12727273 
0.88 



-3.2865187 



-187.64977 

949.29028 

-81.427874 



1.19 

864.31252 

-4439.6671 

321.87224 

-8402.0232 

-18.761161 

-.071698139 

0.99 
1.38 

63.189386 
1064.1109 
-132.91606 

1.00 

8.9650248 
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CONSTANT 
CONSTANT 
CONSTANT 
CONSTANT 
CONSTANT 
CONSTANT 
CONSTANT 
CONSTANT 

CONSTANT 
CONSTANT 
CONSTANT 



KM3 
KM4 
KM5 
KM6 
KM7 
KM8 
KH9 
KMHAC 

KAl 
KA2 
KAHAC 



-10.862163 
346 .68334 
-3818.3687 
54.348067 
-69.482520 
7878 .3217 
-?032.7e86 
4584.2952 
-20049.737 
10.M.4857 
-3972.3512 
20.761223 
-29 .895752 
2. 3778853 
-0.69989459 
U. J9 
n.yj9 
1.06 
0 007 

0.0025286890 
0.006724987 
-.00012021482 
090 



5.2689872 
-2628.0014 . 
1152.1293 . 
-175.99694 , 
0.0 

-25448.989 r 
0.0 

-18887.324 , 
5162.3900 
-6837.4301 f 
1395.9822 
-56.326B84 r 
5.4463004 
-2.2737848 r 
0.090931133 
0.81 

0.946 # 
1.50 

- .0026504608,' 
-.0011479416 
-.0024994776, 
-.00000806355 



4750.3697 

192.95224 

27401.966 

29191.895 

9583.8249 

61.499434 

1.8901715 

0.86 
1.02 

-.00090103102 
.00071838136 



LOGICAL 

CONSTANT 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

ARRAY 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

INTEGER 



INTERA, COSTLM 
INTERA = . FALSE . 



COSTLM 



XNOW( 1600) , 
VXNOM( 1600), 
RHONOM( 16O0) 
RNGNOMt 1600) 
AZMNOM( 16O0 ) 
AZDNOMf 1600) 
VWXPN( 1600} , 
RGMEAS( 16O0) 
A2MEAS( 1600) 
HMFLAG( 1600) r ^ 
1600*^0 



XNOM 
VXNOM 
RHONOM 
RNGNOM 
AZMNOM 
AZDNOM 
RGMEAS 
AZHEAS 
HMFLAG 
COSTTH 
COSTI 
LVX « 



50 



1600*0. 
1600*0. 
1600*0. 
1600*0. 
1600*0. 
1600*0. 
1600*0. 
1600*0. 

0. 9, 

1 . E20r 

. LVY « 



YNOM( 1600) 
VYNOM( 1600> 
VWXNOMr 1600 
RRTNOMI 1600 
ELVNOM( 1600 
ELDNOH( 1600 
VWYPN( 1600) 
RRMBAS( 1600 
ELMEAS( 1600 
MEASW( 1600) 
, YNOM 
, VYNOM 
, VWXNOM 
, RRTNOM 
, ELVNOM 
, ELDNOM 
, RRMEAS 
, ELMEAS 
, MEASW 
PTOL • 
TFINAI* • 
50, LV2 = 



. FALSE . 
ZNOM( 1600) 
V2N0M( 1600) 
VWYNOM( 1600) 
TIMNOM(1600) 



1600*0. , 
1600*0. , 



ZNOM 
V2NOM 



1600*0., WTYNOM 



- 1600*0., 
« 1600*0. 

- 1600*0. 
« 1600*0. 
= 1600*0. 
= 1600*0. 
0.1, BTOL 
90.0 

50 



TIMNOM 



« 0.9 



1600*0. 
1600*0, 
1600*0. 
1600*0. 



NWEIGH, NMEAS 
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CONSTANT 
CONSTANT 

ARRAY 

INTEGER 

INTEGER 

INTEGER 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 



CONSTANT 

ARHAY 
CONSTANT 



CONSTANT 
CONSTANT 



NWEIGH ^ 
NMEAS = 

COSTVX C 4 
NCOST(40 
ILAYER 
RUNFLG 
RUNFLG - 
COSTVX = 
NCOST = 
AX.TLAY ^ 



360, MXWDSL « 0.02, DBFl-AG = 1- 
230 

0), C0STVY(40)r COSTVZ(40), AI,TIAy(20) 
STIAY, PSTIAYr KOONT, HAXIAY. NLAVER 



1, NLAYER - 1. lAYER « 1000, 
40*0, COSTVY « 40*0, C0STV2 
40*0 



TCHNGO = 

RHOI,IM( ? 
RHOLIM = 
1.4245, 
0.90 r 
0.45 - 
0.45 . 
1.0243- 
0.60 
0.25 r 
0.00 , 
VWXLIM = 
VWYIiIM « 



0.0 r 1.0E3, 2 

5.0E3, 6-0E3r 7 

J.0E4r l-iE**^ i 

i.5E4, 1.6E4r 1 
0.5 



.0E3r 3.0E3r 

.0E3r 8.0E3, 

.2E4, i.3E4, 

.7E4r 1.8E4r 



40*0 

.0E3r 

.0E3, 
.4E4, 
.9E4 



0,2), VWXLIM(20,2), VWYI.IM( 20 , 2 ) 



1.40 
n . 80 
0.45 

0 .45 

1 .00 
0.60 
0.00 
0.00 

20* 
20* 



20, 
20 r 



1-20 
0.80 
0.45 
0.45 
0.90 
0.50 
0.00 
0.00 
20*-- 
20*- 



20 
20 



1.00 

0*60 

0.45 

0.45 

0.70 

0.35. 

0.00 

0.00 



0.90 
0.50 



45 
0.45 
0.60 
0.30 
0.00 
0.00 



IF fT.EQ.0-0) ALT = ALTIC 

^TLAY « INT( AI,T/LAYER ) 

OLDALT = ALT 

MAXJUAY - 5TLAY 

JLAYER « STLAY 

CONNT = COUNT - 1 

KOUNT » 0 



+ 1 



IF (T-NE.0.0) GO TO 99997 
Cons^anrs " 



A80PI 
A160PI 



m 8.0 / PI 

16.0 / PI 



- Set seed number " 

INTEGER RNSEED. ^SEEO 
CONSTANT RNSEED « 2234567 
GSEED « RNSEED 
GAUSI(GSEED) 
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Plot initialiratiar " 
COUNT « 0 

Descent flag at T=0 is false 
DSCNTO « .FALSE. 

Initial conditions for integrators 

" usina B?IL FCT Data "^^raiectory 2 : 
Charoe 08 r vmoTtle « 690.8 
Elevation = 53.44 deg*950 MILS, SPIN 



220 REV/SEC=79200 DPS ' 
300 REV/SEC«108000 DPS' 



CONSTANT ELVDIC 
CONSTANT ELVIC 
CONSTANT UIC 
CONSTANT XEIC 



n.O , SPIN 
53.44, ROLL 
690. 8r VIC 

0.0 r YEIC 



79200. r YAWD •0.0 

0.0 , YAW - 0-9 

0-0 , WIC * 2-2 

0.0 , ALTIC • 0*0 



TlIEIC 
PHIIC 
PSIIC 



-ELVIC / 0*^0 RAO 
ROLL / D^roRAD 
YAW / OTORAD 



THDTIC « -ELVDIC / OTORAD 
PHDTIC ' SPIN / OTORAD 
PSDTIC « YAWD / DTORAD 

Proiectile characteristics 



CONSTANT 
CONSTANT 
CONSTANT 
CONSTANT 



MASS 
IX 
lY 
L 



43 .09124 
0. 1477 
1.8009 
. 155 



PI*L*L*0.25 



Mass ( 95 lb ) 

Inertia about x-axis 

Inertia about X-axis 
Diameter 

Section area 



X kg 

5 kg-ni**2 

I kg-ffi**2 

: m 

: in**2 



Constants 



TWOPI - 
L2 

AOHASS- 
IXoIY = 
ALoIX * 
ALOIY « 
IXIX * 
IY2 
IY4 

ALFMAX* 



2. ♦ PI 
L ♦ 0.5 
A / MASS 
IX / lY 
A ♦ L 
A ♦ L 
IX*IX 
2 . ♦lY 
4.*IY 
0.0 



/ IX 
/ lY 



Initialize wind velocity and angle of attack 
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VWBBX 
VWEEY 

SINTHE 
COSTHE 
SIHPSI 
COSPSI 

TJAEB 
VAEB 
WAEB 

l»W 
vw 

WW 
VTW 

ALFA 
EPS 
BPSSQ 
VWoVTW 
WWoVTW 
RHO 
QBAR 
TEMPK 
VSONO 
HACH 

UPEE 
VPEE 
WPEE 

A2IM 
ELEV 
XLOSE 
YLOSE 
ZhOSZ 
RRT 
RNG 

DV2F1 
DVXFl 
DVYFl 



VWXTABf ALTIC) 
VWyTAB( ALTIC) 

SINfTHFtC) 
COS(THEXC) 
SIN(F?IXC) 
COS(PSIIC) 

VWEEX * COSTHE*COSPSI 
-VWBEX * SINPSI^^^^^^^ 
VWEEX * SINTHE ♦COSPSI 



VWEEY * COSTHE*SINPSI 

VWEEY * COSPSI 

VWEEY * SINTHE*SINPSI 



UIC - UAEB 
VIC - v;vEB 
WIC - WAEB 

SQRT( tW*UW+VW*VW+WW*WW) 

ACQS (tw/ VTW) 
SIN (ALFA) 
EPS*EPS 
VW/VT%* 
WW/ VTW 

KRHO ♦ RHOTAB(ALTIC) 
0 . 5 *RHO*VTW*VTW 
KTEMP * TMPTAB{AI»TIC) 
C3 ♦ SQRT( TEMPK ) 
VTW / VSOND 



-UIC*SINTHE * uwoxrx*. 

PSIIC 
-TKEXC 
COS ( ELEV ) ♦COS ( A2 1 K ) 
COS{ELEV)*SrN(AZIM) 

UPEfi^XLOSE ^ VPEE^YLOSE + WPEE^ZLOSE 
0.0 

0.0 
0.0 
0.0 

^^IeTFLgS J*"hr?5utlne RADARSET will allow the user no modify 

to 0 and xxRDAR to nominal vaiues. The feedback 
liir?4ut!SS R^ASSfili'^wSii set HITMSS based on the 



. RADFLG-0! 



RADFLG«ls 
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RADrLG»2 1 



KSIGMA: 



radar equation, xxNOIS to 0 And '^^^J^^^^^^S'^i?^^ « 

The feedback signals are DX, DY and DZ. 
^hi routine RADARMEAS will set HITMSS based on the , 
nSi" "»^rMhold KSIGMA. The feedback signals are ^ 

XERRr ^'^^^^ ^^'^ 2ERR. 
noisi threshold (# of sigmas) 



for hit/mias decision' 



INTEGER SETFLG 
CONSTANT SETFLG 
CONSTANT KSIGMA 



1, MITMSS • l.f RADFLG • L 



IF (RONFLG.BQ-1) CA.L3U RADARSET( RWSEBO, SETFLG) 



PGNOIS 
RGRDAR 
RCSDB 

ADMS 
EDMS 



0.0 
0.0 
0.0 



RRNOIS « 
RRRD^R = 
STND© = 



ELNOIS 
ELROAR 



0*0 
0.0 



$ AZNOIS 
$ AZRDAR 



0.0 
0.0 



0.0 
0.0 



Backward rectangular integration of range rate 

INTGRR « 0.0 



INTEGER NAPEX 
NAPEX = 5555 

iNOFF-w-iiF, ignored 
on. output page 46; 
END S - of INITIAL 



on inuut line 648 of page 27 of file "ANALYSIS : (ARMAL 



DERIVATIVE SLOW 



Simulation control parameters 

NSTEPS ?f7?J " k 

ALGORITHM lALGl = 3 , 

SaXTERVAL MAXTl =0.01 

MINTERVAL MlNTl = 0.01 

Simulation termination conditions 

LOGICAL OVTLOW 

Overflow prediction 
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nvFt^OW . (ABS(U0OT).<5T.0DOT0O) .OR. { ABS{VOOT) .GT. VOOT99 ) 

Atmospheric model, mci inpu« for MET» closed loop e.rmation" 
PROCEDURAL {DRHOFB, nvwxFB, DVWYFB - ) 
IF ( .NOT.CLMET) GO TO NOCI. 

CONSTANT KXER -05, KVXBR - 10, KAXER - 0. S "VHX estim loop" 
OVWXFB «-KXBR*DSF - KVXER*DVXF 

CONSTANT KTER -.05, KVYER ^ 10, KAYER - 0. S "VHT est:i« loop' 
DVWYFB »-KYER*DT^ • KVYER uvxr 

NOCI/. -CONTINUE 
END $ -OF PROCEDURAL - 
CONSTANT KRHO - KVWX . 1., KVWY - 1.. KTEMP - 1. 
CONSTANT VWEE2 • 0- 

•^"^ SiS?Ci;;E?"-?RS^Ft^^o!o, : RHOCORC AI.T , 

^^•^^ " ^^^"'^O^l^Ft^'^lo, : <rWXCOR( A1.T , 
^^^^ " ^CtAEr^W^FS^olo, : ^COR( Al-T ) 
TEMPK - KTEMP * TMPTAB( ALT ) 
. Aerodynamic coefficients 
PROCEDURAL 

speed of eound and Mach number 
CONSTANT C3 « 20. 04*8 

VSOND - C3 ♦ SQRT ( TEMPK ) 
MACH - VTW / VSOND 

Spin deceleration coefficient 
FKA( KA, MACH, KAMAC, KAl, KA2) 
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CLP = - KA • A160FI * CLPFAC 
•- Pitching moment coefficient 

rKM( KM, MACH, KMMAC. KMl , KM2, ICK3, KM4 , KM5 , KM6 , KM7 , KMB,KM9, 
CMA « KM • A80PT 

Damping in pitch momenl: coefficient 
CONSTANT CMQ « 0.0 

Magnus mowent co-tficient 
CONSTANT CNPA =0.0 

Axial force coefficient 

CXO » KDO * A80PI ; CXOFAC 
CX2 = KDA ♦ A80P1 ^ CX2FAC 
CX = CXO + CX2 * EPSSQ 

Normal force coefficient 
CONSTANT KIA «= 5.0 

FKLOC KLO. MACH. KLOMAC, KLOl, KL02, KL03, KL04 , KL05) 
NO = KLO * A80PI 

: ^^0*/'S!^ ePSSQ ) * LFTFAC 
Magnus force coefficient 
CONSTANT KF » O. 13 
CYPA « KF * A160PI 
END $ " of PROCEDURAL " 
- Precession/nutation frequencies 

MU « QBAR*A*L*CMA 

DISC = IXIX*P2 - IY4*MU 

SQdisc - SQRT(ABS(DISC)) 
PROCEDURAL 

IF (DISC.LT.O.^ GO TO LLl 

OMEGA! « (IXP * SQdisc) / IY2 
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0MEGA2 « (IXP - SQdiSC) / IY2 

GO TO LL2 
LLX. .CONTINUE 

OHEGAl • 0.0 

OMEGAi « 0.0 
LL2. .CONTINUE 

END S ■ of PROCEDURAL" 

Gyroscopic stability factor " 

END " " of PROCEDURAL- 
Filtering velocity evror for closed loop MET estimates - 

CONSTANT TAUX - 0-1. TAOY • 0.1. TAUZ - 0.1 

DXri » RSW( RADFLG.EQ.2 , XERR • ) 
DYFl - RSWf RADrLG.EQ.2 , YERR , DT ) 
DZFl - RSwl RADFLG.EQ.a , ZERR . OZ ) 

DXF = REALPLf TAOX . DXFl . 0- ) 
DYF - rIaLPlI taut . DYFl . g. ) 
DZF - REALPLC TAOZ . DZFl , 0. ) 

CONSTANT TAOVX - 2-. TAOVY - 2., TAOVZ - 1. 

DVXFl = RSW( RADFLG.|Q.2 , ^XERR , DVX ) 
nVYFl » RSWI RADFLG.EQ.2 , ^YBRR , Dvy ) 
PVZFI - RSW( RADFLG.Eg.2 , VZERR , DVZ ) 

DVXF - REALPL( TAUVX , DVXFl ' g- } 
DVYF » REA1.PL( TAUV> , gj^^} ' n { 
OVZF - REA1.PL{ TAUVr. , DVZFl , 0. ) 

*'°^^-o;'J«tpGt ilSnt on inout line 818 of page 27 of file "ANALYSIS : fARMAI,. 
END S - of DfeRfVATIVE SLOW - 



OBRIVATIVB^FAST^ 
Simulation control parameters 
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NSTEPS NSTF2 - 1 

ALGORITHM ^^SS ' o 

MAXTERVAL MAXT2 : 0.01 

MINTERVAL MINT2 =0.01 

■ Dynamic pressure (N/Tn''*2) 

OBAR = 0.5 • RHO ♦ VTW ♦ VTW 

> Angle of artack (degi and sine of angle of attack 

ALFA = ACOS(UW/VTW) 
EPS ^ SIN < ALFA* 

EP5SQ » EPS*EPS 
vrwoVTW » VW/VTW 
WWoVTW - WW/VTW 

. Aerodynamic force fv^r unit of mass) and gravity (m/8ec-*2) 

QAM = QBAR : ^'•'^MASS 

PL2V « P * L2 / VTW 

FbxDRG - -QAM^CX 

FbxG = G*SINTH^. 

Fbx « FbxDRG * FbxG 

FbyNF = -QAM^CNA^VWoVTW 
FbyMAG = QAM*PL2V;cypA*WW0VTW 
Fby « FbyNF * FbyMAG 

FbzNF « -OAM*CHA*WWoVTW 

FbzMAG = .QAM*PL2VCYPA*VWoVTW 

FbzG = -G^COSTHE ^ 

Fbr « FbzNF * FbzMAG + FbrG 

Derivative in body axes of projectile linear velocity relative 

to the Earth (in/sec**2) 

finOT « Fbx - THEDOT*w PSIDOT*COSTHE*V 

^OT » Ftoy - PSIDOT * ( COSTHB*U * SINTHE^W ) 

WOT « rbz - PSIDOT*SINTHE*V + THED0T*U 

H Projectile linear velocities relative to the Earth along 

the body x,y,x axes (m/sec) 

U = INTEG ( UDCT,'JIC ) 
V « INTEG VD0T,V1C .) 
W « INTEG ( WOOT.WIC ) 

VT « SORT ( U*U * V*V * W*W ) 
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projectile linear velocities along the Earth axes (m/sec) - 

VPEE » U*C55SJi "-w**-' ^ W»COSTHE 

?RSiE°."sJw(^BEM,rEE * VPEE^VPEE * WPEE*WPEE ) 

Projectile position in the Earth axes (m) 

»F » INTEG ( OPEE. XEIC ) 
YF » INTEG VPEE. YEIC 
ALT - INTEG } WPEF-, ALTIC ) 

Ranges (m) , line of sight unit vector, range rate <m/8) " 

j^ESiE^?-Ji°'/"-n"°" 

YLOSE - YE / RNG 
TLOSE » ALT / RNG 
AZIM - ATAN2f YLOSE, XLOSE) 
ELEV - ASIN(ZLOSE) 
GOTO LOS 2 
LOSl.. CONTINUE 
AZIM - PSIIC 

XLOSE " COS (BLEV)*COS( AZIM) 
YLOSE = C5S BLBvi*SIN(AZIM) 
ZLOSE - SIN(ELEV) 
AZIMD » 0.0 
ELEVO » 0.0 

^^f^i-giS^SIsE * vpee;ylose * wpee-zlose 

IF (RNG.EQ.0.0) GO TO L0S3 

^?^T^OSEL?l4^1^E-4) ELBVD - .WPEE-RRT^ZLOSE) /RNG/COSEL 

ilik°'."t^iB-Rli-^i.isl^^^^ > 

LOS 3. -CONTINUE 
END S - Qf PROCEDURAL 

- Aspect angle (rad) 

cosASP - BooNO( .1.,!.. ^j;gii:^8i?!H:i?glii ::: 

- ZLOSE*SINTHE ) 
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ASPECT - ACOS( COSASF ) 

Aerodynamic momenta uer unit of inertia ( N«,n) / {Kg*m-2 ) 

- "i^*^^ instable p.v-iectile (needs spin to avoid tumbling) : 

- -1 : Stable ptd.i«":tiie 

CONSTANT MSIGN = 1- 

OALIY » QBAR * M.cli 
0ALiy2 = QALiy*L2 

Mtx = qalix*pl:v*clp 

MbvSM » MSIGN*QALIY*CMA*WWoVTW 
MhvDM = OALIY2*'-MQ*THEDOT/VTW 

Projectile angular ».ccelerations (Euler angles) (rad/sec**2) 
PHIOD - Mbx + P?IOO*SINTHE + THB0OT«PSIDOT»C0STHE 

miZ' : S?'5°S;5ioI"rKgSIcOSTHE.FSIDOT-Sl».THE 

pHdD^ ' 'P*gg;'':'';HE6OT*2.*PSIDOT*Sl«THE * Mgyro2 ) / COSTHE 

11^01 I f ^^^.S£^a^??I«^;sioS;^i?^i^HE :^j?s;!^rf /'costhe. 

projectile angular velocities (Euler angles) (rad/sec) 

THEDOT = IKTEG { THEDO, THDTIC ) 
PHIOOT - l»r/EG PHIDD, PHOTIC ) 
PSIDOT - INTEG ( PSIDD, PSDTIC ) 

Projectile attitude angles (Euler angles) (rad) " 

THE - INTEG ( TttEOOT. THEIC ) 

PHI - INTEG PHTDOT, PHIIC ) 

PSI - INTEG ( PSIDOT, PSIIC ) 

PROCEDURAI, ( - THE, VHI, PSI ) 
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Ttr I TMF tST. FT > THE - THE - PI 

\l \ Si -px i the - the + pi 

, OUT PI \ PHI « PHI - PI 

! Ihi :S: Ji 1 phi - - pi 

rr # P^T GT ^\ ) PSI " PSi - " 

'i^ i 111 :S: -PI i PSI - PSI * PI 

END S •■ of PROCEDURAL " 
Sine /cosine o£ Euler angles 

SINTHE • SIN ( THE ) 
COSTHE « COS ( THE ) 
SIMPSl « SIM ( PSI 
COSPSl « COS ( PSI > 

projectile spin rate about the longitudal axis (rad/.ec) " 

P » PHIDOT - PSIDOT*SIHTHB 

. "„d velocity relative to the Earth in body axes (m/sec) " 

UAEB . VWEEX * COSTHE-COSPSI ^ VW||T ^ CO|THE*SINPSI 

^^E^ : '"^^ : ifS?li*cosPsi : ^iv • i???lE*siNPsi 

.. proiectile velocity velative to the atHiosphere in body axes : 

(m/sec) 

uw = U - UABB 
VW « V - VAEB 
WW s W « WAEB 

.. Total projectile velocity relative to the atmosphere (m/s) " 

VTW - SORT ( uw*uw * VWVW * WWWW ) 
.. Projectile hitting ground 

PROCEDURAX. 



LOGICAL COHD4 
CONSTAMT VZLOW « 10 0 

C0N04 - (WPEB.LE.VZLOW) .AND. ( RUNFLG . EQ . 4 ) 
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TEFMT( ALTCND .OR. C0HO4 ) 
END S "OF PROCEOORAL" 

""^'^-'on'lliv^^^ plSe'^tt on i^.-.t line 1041 of page 27 o£ file "ANALYSIS :( ARMAL 
END S" of oIrIVATIVE FAST - 



Layer index 

FSTiIy - INT( OLDALT/ LAYER ) + 1 
MAXLAY - mLo{ MAXLAY , ILAYER ) 
Descent " 

LOGICAL DESCNT , CHANGE , DSCNTO 

. Simulation termination conditions 

LOGICAL ENDTIM, ENOLAY, ENDFLG, EMOCLO 

^S^l?^? 515199 : I'lil :^S?l9 : tho , WD0T99 - 1.E20 

Max Sim time reached 
ENDTIM « (T.GT.TSTP) 

End of altitude layer reached 
ENDLAY = (IIAYER.GE.CSTLAY^NIAYER)) -AND. ( RUNFLG . EQ • 2 ) 

End of closed loop estimation 
ENDCLO . (RONFLG.EQ.4>.A«D.(COUHT.GE.NAPEX-1.0R.COUHT.GE.NMEAS-1> 

ENDFLG - EHDTIM .OR. DESCNT -OR. ENDLAY .OR. OVFLOW .OR. ENDCLO 
TERMT( ENDFLG ) 
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SAVE D^TA FOR PLOTTING USING PLOTIOI' 



INTEGER 
CONSTANT 

Ar.rMdg « 
EliVdeg = 
PSIdeg « 
PHIdeg « 
ALFdeg * 
SPJNdg = 
FREUl = 
FREQ2 = 
ALFMAX - 
THEdeg « 
ASPdeg = 



DATLEN , COUNT 
DATLEN « 1600 



LASTC 
DUMHY 



A2IM ♦ 
ELEV * 
PSI * 
PHI * 
ALFA * 
P 

OMEGAi 
OMBGA2 



PTOPAO 
DTtlRAP 

DT'.'RAD 
PTHRAD 
OT'.JRAD 
Twnpi 
Twnpi 



AMAX1( ALFMAX , ABS( ALFdeg) ) 



THE 
ASPECT 



pTORAO 
* PTORAD 



' Smooth wind computation 



''OMSTANT 
CONSTANT 



VWXAl 
VWYAl 



ft rruvKn - n VWXHl « 5000, VWXHO » 1^95 
8: WYAO = S: ^WHl - ioOO: VWTHO - 1000 



ShX » (ALT - VWXHO J /(VWXHl 



VWXHO) 



ixh^^lir5«§*SL°;g^iSii^2i^,Ai*6HX*DHX*DHX 

SSI- : T^^^-V^S?/(VWyHl - VWYHO, 

IF (ALT-LT-VWYHO) VWVp = VWYAO 
IF ALT-GT-VWYHl) WYp = yWAl 



COUNT = COUNT + 1 

IF (WPEE.LE.VZLOW .AND. RUNFLG.EQ.l 
NAPEX « COUNT- 5 



.AND. NAPEX. EQ-5555) 



IF (RUNFLG.NE.l) GOTO NOSAVE 
TIMNOM( COUNT) = 



XNOM( COUNT) 
yNOM( COUNT) 



= KZ 
= iF 
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ZNOM( COUNT) 



ALT 



VXNOM( COUNT) * U^|f 
VYNOM{ COUNT) = '^IH 
VZNOM( COUNT) = WPEE 

RHONOM( COUNT) = ^^2r^„ 
VWXNOM{ COUNT) = VWEEX 
VWYNOM( COUNT) = VWEEY 

RNGNOM( COUNT) = «NG 

RRTNOM(COUNT) - ^PT 

AZMNOM( COUNT) = A?iw 

ELVNOM( COUNT) = ^•Vfiir^ 

AZDNOMf COUNT) = ^J l^^ 

ELDNOM( COUNT) = CL-EVD 

CALL RADAWlEAS(PAnFLG. 

RGRDARr 

RGNOISr 

LOGICAL NOHM S CONSTANT NOHM 
IF (NOHM) HITMSS = 1- 

CONSTANT ADSIGO = 0.. ^SI^gI - 0* 

igSlIf V?!'""" ^D5iGl*RNGf "GiuSS(0 
IdSoiI - - EDSIG1*RNG) * GAUSS(0 

I F ( COUNT . NE . 1) 00 TO RADAR 



ASPECT, 
KSIGMAr 
RRT, 
RRRDARr 
RHHOI5 , 



HITMSS r RCSDB, 



ELEV, 
ELRDARi 
ELMOIS I 

. FALSE 



AZIM, 

AZRDAR, 

AZNOIS) 



$ "?<SIG<?' 

$ •■?<SIG<?' 
.0, ADSIGO) 
.O^EDSIGO) 



HITMSS 
RGRDAR 
RRRDAR 
ELRDAR 
AZRDAR 
RGNOIS 
RRNOIS 
ELNOIS 
AZNOIS 



1 

RNG 

RRT 

ELEV 

AZIM 

0.0 

0.0 

0.0 

0.0 



RADAR. .CONTINUE 
HMFLAG( COUNT) » HITMSS 



1:. HIT, 0-OOls MISS" 



IF ( HITMSS. EQ. .r- ) 
RGMEAS( COUNT) 
RRMEAS( COUNT) 
AZMEAS( COUNT) 
ELMEASf COUNT) 
GO TO HIT2 

HITl. .CONTINUE 



GO TO HITl 
^ RNG 
= RRT 
« AZIM 
= ELEV 
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RGMEAS( COUNT) 
RRMEAS{ COUNT) 
A3MEAS( COUNT) 
ELMEASf COUNTS 
HIT2. .CONTINUE 



RGRDAR 
RRRDAR 
AZRDAR 
ELRDAR 



IF (COUNT.GT.l) TMTGRR = JNTGRR^-^g . 5^CI^ 



VWXPN (COUNT) = VWXp 
VWYPN (COUNT) = VW\p 

NOSAVE. •CONTINUE 

Cost function computation 

IF (CHANGE) KOUNT » " 

KOUNT » KOUNT * 1 

XNHMO » XNOM(CO"NT) 
YNOMO - YNOMCCOHNT) 
SNOMO « ZNOM( COUNT) 

VXNOttO « VXNOM( COUNT) 
VYNOKO = VyNOM(COffNT) 
VZNOMO = V2NOM( COUNT) 



DX 
DY 
DZ 

DVX 
DVY 
DVZ 

RGMS 
RRMS 
AZMS 
EIiMS 



XE - XNOMO 
YE - YNOMO 
ALT - ZNOMO 

OPEE - VXNOMO 
VP EE - VYNOMO 
WPEE ~ VZNOMO 

RGMEAS( COUNT) 
RRMEAS( COUNT) 
AZMBAS( COUNT) 
ELMEA5( COUNT) 



INTEGER NANGLE S CONSTANT NANGLE - 5 
IEm^'^S^ ( liMslAtSiki 1 c60irr.NANGLE ) ) /CINT/NANGLE 
i^M^^^'^lif^slffiilc^^^^^ ) /CINT/NANGLE 



EDMS 

5INAZ 
C05AZ 
SINEL 
COSEL 



SIN(A2MS) 
COS (AZMS) 
SIN (ELMS) 
COS ( ELMS ) 
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XMEAS 
YMEAS 
2MEAS 

VXHEAS 

VYMEAS 

VZMEAS 

XERR 
YERR 
ZERR 

VXERR 
VYERR 
V2ERR 

RGEPR 
RRERR 



RGMS 
RGMS 
RGMS 



♦ COSEL * 

♦ COSEL * 

♦ SINEIi 



C0SA2 
SINAZ 



- RRHs*cosEL*siNA2 - iiHi:sg«i:i5gg:iiSI 

= RRMS*STNEL * RGMS ♦EDMS* COSEL 

= XE - XMEAS 
= YE - VMEAS 
= ALT " ^MEAS 

= UPEE - VXMEAS 
= VP EE - VYMEAS 
« WPEE - VZMEAS 

a RMG - PGRDAR 
= RRT - RRROAR 



IF (RUNFLG.EQ. 1 » GOTO NOCOST 

XC0ST2 
YCOST2 
2COST2 

TIC0ST2 
VC0ST2 
WC0ST2 

LCOSTX 
LCOSTY 
LC0ST2 

COSTVX(ICOST) 
C05TVY(1C0ST) 
COSTVZ(ICOST) 

IF (OVFLOW) COSTVX(ICOST) = 12|22}" 

IF OVFLOW COSTVY(ICOST) = 

IF OVFLOW COSTVZ(ICOST) = 998001. 



= XERR*XERR 
= yERR*YERR 
« ZERR*ZERR 

= nvX^DVX 
e nVY*DVY 
= DVZ»OVZ 

= XC0ST2 * LVX*UC0ST2 
= VC0ST2 + LVY*VC0ST2 
= ZC0ST2 + LVZ*WC0ST2 

( ( KOUNT- 1 ) *C0STVX ( ICOST ) 
KOUNT- 1 ) ♦COSTVY ( ICOST ) 
KOUNT- 1 ) ♦COSTVZ ( ICOST ) 



LCOSTX 
LCOSTY 
LC0ST2 



) / KOUNT 
) / KOUNT 
) / KOUNT 



COSTX 
COSTY 
C0ST2 

NCOST( ICOST) 



- SQRT(COSTVX( ICOST) ) 
= SQRTCOSTVY (ICOST ) ) 
= SQRT ( C0STV2 ( ICOST ) ) 

= KOUNT 



Feedback loops fov control/identification of METs parameters 
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CONSTANT CLMET = -FIVLSE. 



IF (.NOT. CLMBT) GOT'"' MQCOST 

ARRAY AI.TC(160CM PMOC(1600, VWXCflSOO), VWTC(1600, 

CONSTANT J^TC = \fOo:Oy RHOC _ 1|00^0. 



roNSTANT VWXC = •'^00*0. 
ALTC( COUNT) = AUT 
rSoC COUNT - OR"OF| 
VWXC (COUNT) = pyw$n 
vwyc( COUNT) = DVW-VFB 

NOCOST. .CONTINUE 



irr ALT .LT. 0.0 » «3OT0 NOSAVl 

RHOCO = RH0C0R(A1.T) 
VWXCO » VWXCOR(ALT) 

VHYCO « VWYCORCALT) -,«.,c n O i 

2H0I.E = RSW( HlTHSS-EQ.l. , ZMEAS , 0.0 ) 

CALL SAVEDATAS ( Pf^J^EN, COUNT ,1 , „ ^ BLVdeg, I 

PSXdegl PHIdegl XE , YE , ^ . r • 

ALFMeg QBAR , VTW , XMEAS , YMEAS , • 

ZMEAS , VWXCO , VWYCO , COSTX , UW ) 

CALL SAVEDATAS ( OATLEN, COUNT , |1 ; , THBdeg, 

^^Mrr : SQRT(2C0ST2) SQRT(WC0ST2) , • 

COSTZ , ELN01S*DT0RA0, ..©pf 

VT ' UAEB . VAEB , WAEB , UPEE 



VPEE 



WPEE , Fbx , Fby , Tbz 



CALL SAVEDATAS (DATLEN^COWT , 41^ . ^ ^E, 

XS^rl 'VWEEXT VWXHOM( COUNT)-', DVWXFB, 
?2n2 ; roriNrt YNOMf COUNT ) , 2NOM ( COUNT ) , 

S?cffi) , ^nKcounI; J ™ 

BHONOM( COUNT), CX , EI.MEAS ( COUNT ) DTURAU , 
MACH r INTGRR ) 

CALL SAVBDATAg ( ^"^^|SEAi?COUNT) "rR«EAS ( COUNT) , ... 

azmeaslcount)*dtorad, rgnois, • ^™ zerr 
rrSois, a.2nois*dtorad, xerr ' ' 

RSsDBr-TNDB, THEDD , PHIOD , PSIDD , . • - 
RNG . vwxp , RRT , VWYp , AZlMdg ) 
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i-KTT ClkVCDXTAS f D^TLENr COUNT ,81 ' ^„ 

CALL SAVEDATAb [ ^vNOM COUNT) , DX , DY r D2 , 

Xvy ' DVV DVZ TIMNOM(COUNT) , ... 

SORT(yC0ST2) SQRT(VC0ST2) \ DVWYFB , ... 

RHOCO , HTTMSS , VXERR , VYERR , ... ^^^^ 

V2ERR : WXPN( COUNT) , VWyPN{ COUNT) r ASPdeg ) 
NOSAVX. .CONTINUE 

^^''^-^Jn'output plSe^fo; on inp^t line 1358 of page 27 of file "ANALYSIS : [ARMAL 
END S " Of DYNAMIC " 

TITRMINAL 

IF ( .NOT.OVTLOW) GOTO NOOVF 
WRTTE(20r 123) T 

123,.FORMAT(//,l^.J5H ^^^^^^^ , Overflow prediction r ... 
/ ,1X, lUH > T = r F10.3) 

NOOVF. .CONTINUE 



IF (RUNFLG.NE.l) GOT^* T^ENDl 



COUNT 
LASTC 



COUNT + 
COUNT 



TIMNOM( COUNT) 

XNOMf COUNT) 
YNOM( COUNT) 
2NOM(C0UNT) 

VXNOM( COUNT) 
vyNOM( COUNT) 
VZNOMI COUNT) 



= XF 
« YF 
« ALT 

= UPEE 
= VPEE 
= WPEE 



RHONOM( COUNT) = fHO ^ 
VWXNOM( COUNT) = VWEEX 
VWYNOMt COUNT) = VWEEY 



RNGNOM( 
RRTNOM( 
AZMNOM( 
ELVNOM( 
AZDNOM{ 
ELDNOM( 



COUNT ) 
COUNT ) 
COUNT ) 
COUNT ) 
COUNT ) 
COUNT) 



RNG 

PRT 

AZTM 

ELEV 

AZIMD 

ELFVn 



+ RGNOIS 
+ RRNOIS 
+ AZNOIS 
+ ELNOIS 
+ ADNOIS 
+ EDNOIS 
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ymxw (COUNT) = vwxp 

VWYPN (COUNT) - VWVp 

HMFLAG( COUNT) « HITMSS 

HGMEAS( COUNT) 
RRMEAS( COUNT) 
A2MEAS( COUNT) 
ELMEAS( COUNT) 



n?Wf HITMSS.EQ.l. 

RWf HITMSS.BQ.l. 

PSW( HITMSS.EQ-l- 

HSWC HITMSS.EQ-1. 



RGRDAR r RNG ) 

RRRDAR r R?T >, 

AZRDAR r AZIM 

ELRDAR r ELEV ) 



CAIiL SAVEINF0(1 



KIT r YE 

n . O ,0.0 

CALL NAMEOFPtOT 

CALL PLOTROOTINE tsBrirrQ<;MEAS 
IF (.NOT-NOHM) CALL PROCESSHEAS 



ALT 
0.0 



0-0 
0.0 



0.0 
0.0 



LENDl. -CONTINUE 



IF (RUNFLG.NE.3) GOTO LEND2 



XNOMO 
YNOMO 
ZNOMO 

VSNOMO 
VYNOMO 
VZNOMO 

DX 
DY 
DZ 

OVX 
DVY 
OV2 



« XNOM(LASTC) 
« YNOM(LASTC) 
= ZNOM(LASTC) 

« VXNOM(LASTC) 
« VYNOM{LASTC) 
= VZ^'IOMCLASTC) 

« XE • XNOMO 
« YE - YNOMO 
« ALT - «NMMO 

« UPEE - VKNOMO 
» VTEE - VYNOMO 
= WPEE - V7,NOM0 



YE 
LVZ 



ALT 
DX 



DIMPAC, LVXr - • 
DY f 0.0 ) 



DIMPAC - SQRT( DX*DX + OY^OY ) 
CALL SAVEINFO(ll, XE^. 

CALL NAMEOFPLOT 
CALL PLOTROUTINE 

LEND2. .CONTINUE 

IF fRUNFLG.NE-4) GO T^^^l^^l^ypo TAUVZ , KXERr KVXERr 
CALL SAVEINF0(21, I^er: TAUVY. 0.0) 
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CALI» NAMEOFPLOT 
CALL PLOTROUTINE 
LEND3 . . CONTINUE 



END S " of TERMINAL 

ij^OFF-w-IIf, L ignored 
o 

END 



on output page 62; on inpu 
S of PROGRAM " 



t line 1459 of page 27 of file "ANALYSIS : [ARMAL 



SUBROUTINE NAMEOFPLOT 
INCLUDE 'COMFLOT' 



NAME ( 1 ) 
NAME ( 2 ) 
NAME ( 3 ) 
NAME ( ^ ) 
NAME ( 5 ) 
NAME { 6 ) 
NAME ( 7 ) 
NAME ( 8 ) 
NAME ( 9 ) 
NAME( 10) 
NAME( 11) 
NAME( 12) 
NAME (13) 
NAME( 14) 
NAME( 15 ) 
NAME ( 16) 
NAME( 17 ) 
NAME (18) 
NAME( 19) 
NAME( 20) 



NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAME( 
NAM£( 
NAME( 
NAME( 
NAME( 
NAME( 
NAMEf 



21) = 

22) = 

23) « 
24 ) « 

25) = 

26) = 

27) = 

28) = 

29) = 

30) = 

31) « 

32) = 

33) = 

34) = 

35) = 

36) = 

37) = 
38i = 

39) = 

40) = 



'ALT 

'U 
•V 
'W 

' ELEV 
'FSI 
'PHI 
'XE 
'YE 
'TIME 
'ALFA 
' QBAR 
' VTW 
'XMEAS ' 
'YMEAS ' 
' ZMEAS ' 
'VWXCOR' 
' VWYCOR ' 
' COSTVX ' 
'UW 

'VW 

'SPINRATE' 
'RHO 
•THETA ' 
'TEMPK ' 
•ZCOST ' 
'WCOST • 
' COSTVZ ' 
'ELNOIS' 
'VT 
'UAEB 
'VAEB ' 
'WAEB 
'UPEE 
'VPEE ' 
'WPEE 
'Fbx 
'Fby 
'Fbz 
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NAME f 4 1) 
KAME(42] 
NAME M 3 ] 
NAME (44 ] 
NAME(45: 
NAME f 4 6 : 
NAME (47] 
NAME f 4 8 ] 
NAMEf 49: 
NAME (50] 
NAMEfSi] 
NAME f 52; 
fiAi«Ef 531 
NAME ( 54 : 
NAME (55; 
NAME f 56] 
NAME (57; 
NAMEfSO] 
NAME f 5?] 
NAME (601 

NAMF(6i: 
NAME (62; 
NAME(63] 
NAME (64] 
NAME(65] 
NAME (66: 
NAME (67 : 
NAME (681 
NAME (69) 
NAME (70] 
NAME(71] 
NAME (72; 
NAME (73; 
NAME (74 ; 
NAME (75] 
NAME ( 76 : 
NAME (77 ; 
NAMEC78: 
NAME{79: 
NAME (80) 

NAME(83 , 
NAME ( 82 ] 
NAME(83] 
NAME (84] 
NAME (851 
NAME (86] 
NAME (87] 
NAME (88] 
NAM£(89| 



'VXMEAS* 
'VYMEAS* 
'V2MEAS' 
'MbySM • 
' ZDETECTEO ' 
' DRHOFB • 
'VWEEX • 
'VWXNOM' 

' XNOM * 

' YNOM ' 
'2NOM 

'VXNOH • 

'VYNOM ' 

'VZNOM • 
'RHONOM' 

'CX ' 
' EIiMEAS 

'MACH 

' INTGRR ' 

'SG 

' PGMEAS ' 
' RRMEAS ' 
'AZMEAS* 
'RGNOIS' 
'RRNOIS' 
'AZNOIS' 
' XERR 
'YERR 
'ZERR 
^RCSOB ' 
'STNDB • 
'THEOD • 
'PHIDD ' 
'PSIDD * 
•RANG 
'VWXp ' 
'RNGRAT* 

' vwyp 

'A2IM ' 

: 'VWEEY ' 

: 'VWYNOM* 

r 'OX 

: 'DY 

: '02 

I 'OVX 

: 'DVY 

5 ' DVZ 

s 'TIMNOM- 
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NAME (90) = 
NAME (91) = 
NAME (92) = 
NAME (93) - 
NAME (94) = 
NAME (95) = 
NAMr(96) = 
NAME (97) = 
NAME (98) « 
NAME ( 99 i = 
NAME( 100)= 

NAME( 101 )= 
NAME( 102 )= 
NAMEf i03)« 
NAME (? 04 )« 

INFO( 1) = 
INFCM 2) « 
INFO( 3 J = 
lNFO( 4) = 
INFO( 5) - 
TNFO( 6) - 
INFO( 7) = 
INFO( 8) = 
rNFn( 9) « 
INFO (10) » 
INrO( 11 ) = 
INFO(12) « 
INF0(13) « 
INF0(14) ^ 
INFn(15) = 
INF0(16) = 
INF0(17) = 
INFO( 18 ) = 
INFO (19) = 
TNFO(20) = 
INF0(21) « 
INFO(22) = 
INFO(23) = 
INFO(24) = 
INFO(25) = 
1NF0(26) = 
INFO(27) = 
INFO(28) = 
INFO(29) = 
INFO(30) = 

RETURN 



'YCOST ' 
' VCOST • 
' DVWYFB • 
' RHOCOP ' 
'HITMSS' 
'VXERP ' 
'VYEPP ' 
' V2EPR* 
'VWXPN' 
'VWYPN' 
'ASPECT' 

' RGMEAS2 • 
'RRMEAS2 ' 
'ELMEAS2 • 
'A2MEAS2 • 

'XIMPCO' 
' YIMPCU' 
'2;IMPC0• 



'XIMPC1 ' 
' YIMPCl ' 
ZIMPCl ' 
' L'IMPAC ' 
'LVX' 
'LVY' 
'LV2 • 
' DXIMPAC 
'DYIMPAC 

t r 

'K2ER ' 

'KVZER ' 

' TAUV2 ' 

' KXER ' 

* KVXER ' 

•TAUVX ' 
' KYEH 

•KVYER ' 

'TAUVY • 
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SUBROUTINE SET ACSL PARAM 
RUNFLGi f RADFV- 



VWXTABi 
VWXHOi 
VWYAOi 
KXERl 
KAYERi 
KRRERi 



VWYTABi 
VWXAHi 
VWYHli 
KVXEP i 
K2ERi 



RNSEEDi 

LVXi 

VWXHli 

VWYAli 

KAXERl 

KVZERi 



NLAYERi 

LVYi 

VWXAli 

CLMBTi 

KYERi 

KAZERi 



RHOTABi 

LVZi 

VWYHOi 

TSTPi 

KVYERi 

KRGERi 



5 PASS ARGUMENTS TO ZZCOV 



S 2SC0K 

INTEGER 

REAL 

REAL 

REAL 

LOGICAL 

REAL 

REAL 

REAL 

REAL 

PTINFLG 
RADFLG 
RNSEED 
NLAYER 



RtiNFT^Gi RADFLG Jl. RMSEEDir NLAYERi 
RHOTABi U0)7wk^^^ VWYTABi(40) 
LVXir I*VYir LVZi 
TSTPi 

CLMBTi ^ 
KXERir KVXERi. 
KYERi r KVYERi, ^A\ERi 
KZERi, KVZERX, KA«.ERJ. 
KRGERi r KRRERi 

= PUNFLGi 
* RADFLG i 
= RNSEEDi 
» NLAYERi 



CLMFT 

IF f TSTPi 

KXEP 

KVXER = 
KAXER 

i:yer 

KVYER 
KAYER « 
KZER 

KVZER « 
KAZER « 
KRGER 
KRRER « 



CLHETi 
.GE.O-) TSTP 
KXERi 
KVXERi 
KAXERi 
KYERi 
KVYERi 
KAYERi 
KZERi 
KVZERi 
KAZERi 
KRGERi 
KRRERi 



T5TPi 



RHOTAB(i) - RHOTAEi(I) 
VWXTAbU ' VWXT^RKI) 
VWYTAB i « VWYTABi(I) 
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END DO 

LVX = LVXi 
I.VY = LVYi 
LVZ = LVZi 



VWXHO 
VWXAO 
VWXHl 

\rwxAi 

VWYHO 
VWYAO 
VWYHi 
VWYAl 

RETURN 
END 



= VWXHOi 
« VWXAOi 
^ VWXHli 
= VWXAli 

= VWYHOi 
= VWYAO i 
« VWYHli 
» VWYAl i 



1 KXERO r KVXERO , KAXERO 

i; KAYERO r KZERo , vwYTABo ) 

l' KRRERO , RHOTABC r VWXTABO , VWYTABO ) 



C0STV2O 

KYERO 

KAZERO 



NCOSTO 
KVYERo 
KRGERo 



C RETRIEVE VARIABLES FROM "ZOOM 



ZZCOM 

INTEGER MLAYERO, NCOSToMH) ^^nsrVZoMO) 
PEAL C05TVXO(40), COSTV\o(40), C05TVZO(4U) 

REAL 
REAL 
REAL- 
REAL 
REAL 



KXERo, KVXERO, KAXERo 
KYERO, KVYERo, KAYERO 
KZERo, KVZERO. KAZERo 

rSoT^o ( 4 0?"^XTABO (40), VWYTABo (40) 



IF ( RUNFLG . EQ . 1 ) MLAYERc = MAXLAY 



DO I « 1 



40 

COSTVXO(I) 
COSTVYO(I) 
C0STV2o(I) 
NCOSTO (I) 

RHOTAco(i) = 

VWXTABO ( I ) = VWXTAB ( I ) 
VWYTABoll) » VWYTABdl 



COSTVX ( I ) 
COSTVY(I) 
COSTVZ ( I ) 
NC05T (I) 
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END DO 

rxEP.o 

KVXEP.O 

KAXERO 

KYERc 

rVYERO 

KAYERO 

KSEPO 

FVZEPO 

KASEPO 

KPGFPo 

KRREPo 

rrruRN 

EVV 



* KXER 
= KVXER 
= KAXER 
= KYER 
= KVYER 
= KAYER 
- KZER 
= KVZER 
« KA2ER 
= KRGER 
= KRRER 



" SUBROUTINE SAVENOKINAL ( TILE ) 



S SSCOK 

CHARACTER* 10 FILE 

REAL XNQMsav( 16 00) r 
REAL VXNOMsav( 1600) , 
P.EAL RHONOMsavf 1600) , 
PEAL RNGNOMsav( 1600) 
REAL A2MN0M5av( 1600) 
HEAL AZDNOMsavf 1600) 



VNOMsav( 1600) r 
VVNOMsav( 1600) f 
vwXNOMsav( 1600) < 
RRTNUMsav( 1600) , 
ELVNOMsav( 1600) 
ELONOMsav( 1600) 



ZNOMsav( 1600) 
VZNOMsav( 1600) 
VWYNOMsav( 1600) 
TIMNOMsav( 1600) 



ELMEASsav(1600) 



REAL VWXPNsav( 1600 \ . ^VWYPNsav 1600 ) 
PEAL RGMEASsav 1600), RPMEASsav ( 1600 ) 
PEAL AZMEASsav. 1600) 
REAL HMFLAGsav( 1600) 

INTEGER Ir COUNTsav 
INTEGER NAPEXsav 



r Save current ZZCOM arrays 

PRINT *, 'Save nominal values in x.BIN file' 



NAPEXsav 
COUNTsav 

C PRINT *r 



NAPEX 
COUNT 

COUNTsav ' , COUNTsav 



SUBSTITUTE SHEET 



wo 92/19928 - — PCr/US92/03377 



95 



1 



COUNTS a V 



XMOMsav(I) 
YfiOMsav(I) 
7N' »Msav( I ) 

VXN'JMsav( I ) 

vvrnjMsav( I ) 
V2N0Msav( I ) 
P.HOMOMsav( I ) 
\*V»XNOMsav( I ) 
vwyMaMsav(I) 
PNGNOMsav( I ) 
PRTNOMSavj I ) 

A6MNC>Msav( I ) 
CLVNnM5av( I ) 
AZDNUMsav( I) 
ELDN0M5av( I ) 

VWXPNsav(I ) 
VWYPNsav( I ) 



XNCM( X \ 
YNOM( I ) 
ZNOMr X 1 
VXNOM( I ) 
VYNOM ( I ) 
VSHOM( ) 
RHONOM( T ) 
VWXNOM ( I > 
VWYNOM i 1 ) 
PNGNUM( T ) 
PRTNOW { 1 ) 
TIMNOM { T ) 
AZMNOM [ ^ ) 

ELVNO?^( ^ ) 
AJDNOMf T ) 

EliONOM' T ) 

VWVPN(I) 



PGMEASsav( I ) 
PPMCASsavi I ) 
A2MEASsav( I ) 
ELMEASsav( I ) 
HMFLAGsav( I ) 
END DO 



RGMEAS ( T ) 
RRMEAS ( ) 
A2MEAS ( J ) 
EI.MEAS ( I } 
HMFLAG( T 1 



C Restore Z2C0M from ^ IL.E 

CALL ZZSVRS2( 0, FILE) 
r Transfer saved arrays into new ZZCOM 

L,A5'rr = COUNTsav 
NAPEX = NAPEXsav 



DO 



J = i , 
XNOM( T. ) 
YNOM( I ) 
ZNOM( I ) 
VXNOM ( I 
VYNOM( I 
VZNOM ( I 
RHONOM( 
VWXNOM ( 
VWYNOM( 
RNGNOM( 
PRTNOM( 
TIMNOM( 
A2MN0M( 
ELVNOM( 
AZDNOM( 



COUNTS av 

= XNOMsav( T \ 
= YNOMsav( I ) 
= ZNOMsavf T ) 
) = VXNOMsav( 1 ) 
VTNOMsav(I) 
VZNOMSBvf T ) 
RHONOMsav( I ) 
VWXNOMsav( I ) 
VWYNOMsax*( J. ) 
RNGNOMsav( T ) 
RRTNOMsav( I ) 
TIMNOMsav( I ) 
AZMNOMsa--( X ) 
ELVNOMsavf 1 ) 
AZDNOMsavj I \ 
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VWYPN I = VWYPNsavfl) 
RGMEASsavr T ) 
RRMEASsavf I ) 
AZMEASsavf T ) 
ELMEASsavi r > 
HMFLAGsa"f T) 



KGMEAS ( I ) 
P.RMEAS( I) 
AZMEAS( I) 
TLHEAS ( I ) 
HMFLAG(I) 
END 00 



Fill additional elements 



with last nominal values 



DO I 



= HIN ( 1600 , COUNTsa V. n , MIN ( 1600 , COUNTsav^lO ) 



:nom( 1 ) 

YNOMf T ) 
ZNfJMf I ) 
VXNOMf X ) 
VYNOM? I) 
VJNOMr I ) 
PHONOMt T) 
VWXMOMf I ) 
VWYNOM( I ) 
HNGNOM( I ) 
RRTNOM(I) 
TIMNOM( I) 
AZMNOMfl) 
ELVNOMf I ) 



RGMEAS ( I ) 
RRMEAS ( I ) 
A2MEAS ( I ) 
ELMEAS ( 1 ) 
HMFI*AG( I) 
END DO 



WXPN(I) 
VWYPN ( I ) 



XNOMsav ( cvUNTsav ) 
YNOMsav ( cOTiNTsav ) 
ZNOttsav f » *' lUNTsav ) 
VXNOMsa**^ ' Cni JNTsav ) 
VYNOMsav ( COHNTsav ) 
VZNOWsav ( fjf njNTsav ) 
RHONOMsav r cnMNTsav) 
VWXNOMS a V f COUNTS av ) 
VWYNOMsav( COUNTS av) 
RNGNOMS a V ( COUNTS av ) 
RRTNOMsav ( COUNTS av) 
TIMNOMsav ( COUNTS av ) 
AZMNOMsav ( COUNTS av ) 
ELVNQMs a V ( COUNTS av ) 



vwXPNsav ( COUNTsav ) 
vwyPNsav(COUNTsav) 



RGMEASsa V ( COUNTsav ) 
RPMEASsav f COUNTsav ) 
AZMEASsav ( COUNTsav) 
ELMEAS sav ( COUNTsav ) 
HMFLAGsav c COUNTsav ) 



Save 



modified ZZCO** back into FILE 



CALL 2SSVRS2( Ir FILE) 



RETURN 



^ SUBROUTINE SAVECORRECTION ( FILE . NBPMAX , ICORRECT ) 



S SSCOK 
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?NT?rER NBPMAX ICORRECT, NBPsav, NBREAKP0INT5 , I, 
PeIE^rSoCOR^(3200^^^ VWYCORsav(3200) 
REAL FLNEAR r ALTITUDE 
FEAL SM'^OTHEDC 1600) , CO^^^i^) 
RiXt ivERAGED(1600), WEr^irCO ( 100 ) 
INTEGER NWElGHT/50/ ^ ^^^^ 

DATA WEIGHTU / 25 * 1 . , 25 * " . 5 , 50*0 
INTEGER WEIGHT POWER/0/, IER_LSQ 
REAL RHO_WEIGHTS( 1600) 



.25 / 



r — Tranc:f«r xxxTAB xxxC (with xxx ^ RHO, VWX, VWY) so ^ihat the 

isrim«ed pSra^eter xxx is^in xxxCOR which then can be smoothed 



IF (I CORRECT . EQ - 0 ) THE" 
DO I = 1 r COUNT 



ALTITUDE 
RHOC( I ) 
WXC( I ) 
VWYC( I ) 
END DO ! I 
END IF 



ALTC { T ) 
RHOC ( I > 
VWXC ( I } 
VWYC(I) 



FLNEAR (20,20, RHOTAB , ALTITUDE ) 
FLNEAR 2 0 , 2 0 r VWXTAB , ALTITUDE ) 
FLNEAR 20 20, VWYTAB , ALTITUDE ) 



Sum-up correction values at the different breakpoints 



FLNEAR ( NBPMAX , NBP , RHOCOR , ALTITUDE ) 
FLNEAR NBPMAX , NBP , VWXCOR , ALTI TUDE 
FLNEAR NBPMAX , NBP r VWYCOR . ALTITUDE ) 



IF ( .t CORRECT. NE.O) THEN 
DO J = I , COUNT 

A.LTITUDE = ALTC ( I ) 
RHOC ( I ) =^ RHOC ( I ) 
VWXC ( I ) = VWXC ( I ) ' 
vwyc ( I ) = VWYC ( I ) 
END DO ! I 
END IF 

Smooth-out RHOC with a quadratic polynomial 

xr ( DBFLAG . EO. 1 . ) THEN 

PRINT ' TAR> smooth out RHO ? (l:yes, 0:no) : 

READ(5r*) I 
ELSE 

I « 1 
END IF 

IF (I.EQ.l) THEN 

I = MTN( COUNT , NMEA5 ) 

CALL LS0_POLy2( COEFF , ALTC, RHOC, I) 

iMOOTHEDm'^ COEF-(3. * ^""^l 'Al^C 1 if ^IoIfF ( 1 ) ) 
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END DO i I 

CALI. PLOTSMOOTHEDt SMOOTHED, ALTC, RHOC. COUNT, 

l>R!OT'*;^'*TAR>'^e"p smoothed RHO ? (Ityes, 0:no) : ' 

RBAD(5,*) I 
ELSE 

1 = 1 
END IF 

IF fI.EQ.l) THEN 

^h5c(I) =°SMOOTHeO(I, 
END DO I I 
END IF 

END IF 

r smcoth-out vwx by e^'sraging 

^^?^^''^:^'^*kA>™Sool:h cut VWX ? (l:yes, 0:no) : ' 

READ{5.*) I 
ELSE 

I * 1 
END IF 

IF (I.EU.l) THEN 
200 CONTINUE 

gA£l-"MSU?ESlGE(''A^H^iEO. VWXC, J, WEIGHTO, HWEIGHT, 

END DO 

CALL PLOTSMOOTHEDC AVERAGED, ALTC. VKXC, COUNT) 

^R!^T''*:''°"TAA>''Klep smoothed VWX ? (l:yes, 0:no) : ' 

READ (5,*) I 
ELSE 

1 = 1 
END IF 

IF (I.EQ.l) THEN 
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250 



DO I 1 f COUNT 

VWXC|I) = AVERAGCO(I) 
END DO ! I 

^^^PRTNT * ' TAB> Chgnae weights ? (l:yes, 0:no) : 

HRINT**^* ' ^TAR>'^Num^^^^^ ? NWEIGHT 

PEADfS,*) NWEIGHT 

IF ( NWEIGHT. LE.O) TO 290 

PRTMT * ' TAR> Enter weionts : 

PEAD(5,M WEIGHTO 
GO TO 200 
END IF 



END IF 

2 90 CONTINUE 

C Smcoth-out vwy by averaging 

^S!^;^'-5:^? iAi>lmooth out VWY 7 (l:yes, 0:no) : " 

READ(5,*) I 
ELSE 

1 = 1 
END IF 

IF (I.EQ.l) THEN 

3 00 CONTINUE 

^A£L"Mn;J?^^GE(''Av|RA<^ED: VWYC, J, WEIGHTO , NWEIGHT, 

DO I = NMEAS + 1 , COUNT 

AVERAGED(I) = AVERAGED ( NME AS ) 
END DO 

I limii- max wind variation to MXWDSL (ni/sec)/in 
Cpll tmiTyinVi Xvl^GED, ALTC, COUNT, MXWDSL) 

CAIi PLOTSMO0TKED( AVERAGED, ALTC, VWYC, COUNT) 

i>S!ST''*;^9-TAA>''Klep smoothed VWY 7 (l:yes, .0:no, : 

P.EAD(5,*) I 
ELSE 

1 = 1 
END IF 
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IF r T . EQ . 1 ) THEN 

^;C(I) S°A^RA.iED(I) 
END DO ! I , 
^'--PRINT *, ' TAR> Chenge weights ? (Ityes, 0:no) t 
READ (5,*) ^ „ 
;i.l,^'^.^-^'rlt>'"u.^^^- of weights ? NWEIGHT 

^pn? 250, TwElGHTn(I),I=l,NWEIGHT) 
TO 300 

end' IF 
EHD IF 
330 CONTINUE 

. Load values inte "-^av- arrays 

NBPsav = COUNT 

I = 1 , NBPsav 
' J = l^NBPMAX ! index of independent value 

vwxCORsav(I) » y?^^ ^ 

WYCORsav(I) = yWYC I 
vwYCURsav(J) * ALTL(l) 

values into ar.ays (for extrapolation use, if enough roo. 

I = COUNT * 1 

IF (I.LE.NBPHAX) THEN 

NBPsav = I 
J = I^NBPMAX 



ALTITUDE - 2.*ALTC(C0ttNT) 

PHOCORsav(I) = ^OCOP^av(COUNT) 
P.HQCORsav( J) = ALTITUDE 
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VWXCORsav(I) = VWXC0B5av( COUNT) 
VWXCORsav(J) = ALTITUDE 

VWYCORsav(I) - VWYCORsavf COUNT) 
VWyCORsav(J) = ALTITUDE 

END IF 

r Save correction tables into nominal 22C0M file 

CALL S2SVRS2( 0 r FILE ) ! restore 
NBP = NBPsav 
DO I » 1 , NBP 

J = I+NBPMAX ! index of independent value 

PHOCOP(I) = RH0C0RS9^'( T ) 

RHOCOR(J) = RH0C0Rsav( T ) 

VWXCOR(I) = VWXC0Rsav( ) 

VWXCOR(J) = VWXCORsav(J) 

VWYCOP ( I ) = VWYCORsa V ( I ) 

VWYCaR(J) = VWYCORsav(J) 



END DO ! X 

IF ( ICORRECT. EQ.O ) THEN 
DO I = 1 . 2 0 

RHOTAB{ I ) = 0.0 
VWXTAB ( I ) = 0.0 
VWYTAB(I) = 0.0 
END DO ! I 
END IF 



! erase xxxTAB since their values 
! have been trans fered to xxxC 
! at the beginning 



CALL Z3SVRS2( 1 
RETURN 



FILE ) 



save 



SUBROUTINE PROCESSMEAS 

Thi«^ routine processes f.he radar measurements prior to their use 

to drivo the closed loop estimation algorithm. , ^ « . ^i^«T-c»H 

The rana- ranoe ra^ elevation anale and azimuth angle are considered 
The rang^, "SI time and approximated by polynomials of order P^ORDER 



as functions 
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. a least-s^a.e sense .approximation up .o ..NK.XOH^CXHT, . 
C The purpose IS twofold^ . ^^en HITMSS xs zero, 

^ : 1% iSooth%ur!h2°nSise in the measurements. 



S ZZCOK 



<- IMpLICI''" NONE 
^ REAL. RGMEAS(l) 
r- REAL MEASWf 1) , 
^ TNTEGEB NWEIGH 



i.R«EAS(lW ELMEASd), AZMEASd). TIMNOMd, 



INTEGEP 

P ORDER 
FT 



polynomial order 
polynomial order 



+ 1 



PARAMETER ( P_ORDER.2 , F .L ««.0RDBR*1 ) 



INTEGER 

N DATA 



REAL 



NSEGMENT/2/ . ^ 
I SEGMENT i index 
INDEXlr INDEX2 
I' J 

RNGOPTION/l/ 
NEXTRA 



P COEFfPl) 
TWORKf 1600) 
RGWORK(16 00) 
HRWORK(1600) 
ELWORK(1600) 
AZWORKf 1600 5 
WKARRAY1(1600) 



INTEGER WEIGHT/20/ 



NPOmTS 



number of data , poin« for approximation 

l^l^r El'?SriigSS?s considered in NWEIGH^CI 



indices 
indices 
option index 



Dolynomial coefficients 
work array for txme 
work array for range 
wSrk array for range range 
wSrk airay for eierratxon 
work array for azimuth 
work array 



number of weights for moving a v« ^^^^^^ 
total number of points m sequent 



REAI« 



^iv wPTGHTSfSO) weights for moving average 

DATA AV>E!iHTlHE2i..inM.K5j 



C- 



Farameter 



initialization 
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NWEIGH - MIN( NMEAS , HAFEX ) 

i::::::-£;;r;;rr"sr;iis^sr„F"^'""'""°"""^^^ 

3 polynomial aPPfO*i;r!-i-2! 

100 CONTINUE 

" ^R!S^?;%Uerif r-e segmenrs (min=l) ? default - 
ERim , NSEGMENT 

N5EGMENT « MAX( NSEGwriNT , i , 

END IF 

DO 1 SEGMENT - 1 , NSEGMEVT 

TMnrvl - NWEIGH*(ISEGHENT-1V/NSEGMENT + 1 
-INdIx2 ; {JlJiiiS-ISEGMENT/NSEGMENT 

Make sure there are enough data points 

N_DATA - 0 
PRINT *, 1HDI.X1. INDEX2, NWEIGH 

ir=(^Eil(i,'ES"^') N.DATA . N.DATA . 1 

END DO 

1^ ^,STNf--'''Norinluah data points.'. N.DATA 
IF ( I SEGMENT. EQ.l) THEN - 

i?°f^^Drxi'il^;^Cl4H) THEN 

?r?^EIGS'*ilxO^Nl, THEN , 

i.RINT *, •'""^'t ^^^^l^iSl^ld^Sf^el this plot session.' 

i^^ ?k^aa?"pirlormance to be increased ' 

END 

GO TO 100 
END IF 

^"iNDEXl - INDEX 1 - 1 



C 

200 
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END IF 
GO TO 200 
END IF 

^ orocess ranae measurements 

' PRINT *r ' range' 

DO^^I^- XNDEXl . INDEX? 

i«^^^<^>-^^-^-^TSATA. 1 

SkarrKykn DATA* - 2?iISSS{i| 

TW0RK(N_DATA) * TIMNOM(I) 

END IF 
END DO 

IF fM DATA.LT.3) THEN 

PRINT 'Error, n^data = , n_data 

STOP 
END IF 

CALL LSQ^F0Ly2( P.COEr, TWORK, WKARRAYl, N^DATA) 

DO I ' INDEXl . INDEX2^NE}CTRA 
IT^ { HMFLAG ( I ) . NE . 1 • iJ??^ 
RGWQRK(I) « P LOEF(l) 

RG^ORK(I) - PGiWORK(I) * TIMNOM(I) * P_COEF(J) 
END DO 

^^^RGWORK(I) = RGHEAS(I) 
END IF 
END DO 

RGWORK{I) = P COEFd) 

'"^ ^GSoRKd)'" RGWORK(I) * TIMNOM(I) * P_COEF(J) 
END DO 
END DO 
END IF 

r- Process ranae rate measurements 

PRINT ' range rate' 

N DATA - 0 ,„r.r.v-> 

DO I « INDEXl r INDEX? 

WFARRAY1(N_DATA» » RRMEAS(I) 
END IF 
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jr »N DATA.lt. 3) THEN 

trTnt *r 'Error, n.hata « / n_data 

STOP 

CND IF 

CALL. LS0_F0Ly2( P_COEF, TWQRK, WKARRAYl, N_DATA) 

DO T = TNDEXl , INDEXJ-^MEXTKA 
IF (H«FLAG(I) .NE. 1;^) THEN 
RRWORK(I) = P^COEF(i) 

RRWORK(l!^« RRWORK(I) * TIMNOM(I) + P_COEF(J) 
END DO 
ELSE 

RRWORK(I) = RRMEAS(I) 
END IF 

END DL' 

ir ( T SEGMENT -EQ.NSEGMFNT) THEN 

Df> T = NWEIGH+1 , NAPEX+NEXTRA 
RRWORK(I) = P COEF{i) 

RRWORK(lI^- RRWORKd) * TIMNOM( I ) + P.COEF ( J ) 

END DO 
EHD DO 
END IF 

Process elevat'.'^n angle measurements 

PRINT *r ' elevation' 

N DATA =0 

nn 7 = INDBXl , T.NDEX.^ 

IF (HMFLAG(I) .EQ..I. ) THEN 

N DATA = DATA 1 

WKARRAY1(N DATA) = EEMEAS{I) 
END IF 
END DO 



IF (N DATA.LT.3) THEN 

PRINT 'Error. N^DATA = N_DATA 

STOP 
END IF 

CALL LSQ_P0Ly2( P^COEF , TWORK. WKARRAYl, N.DATA) 

DO I = INOEXl , INC)EX2-NEXTRA 
IF (HMFLAG(I) -NE.l- ) THEN 
ELWORK( I) = F OJEF( 1) 

ELWORK(lI^- ELWORK(I) * TIMNOM(I) P_COEF(J) 
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END DO 

^^^IlWORKCD = ELMEASd) 
END IF 
END DO 

DO I - l^'EIGH-l , WAFEX|NEXTRA 
ELWORK(I) = F COEF(l) 

ELWOWCClT - EL.WORK(I) * TIMNOM( I) * P.COEF(J) 
END DO 
END DO 
END IF 

process arimuth angle measurements 

PRINT ♦r • arimuth 

INDEXl . ^^^f^^r, 

Sk^I^YKN^DATAi = AIMEAS(I) 
END IF 
END D'J 

^,^T£Jio^'^.OATA . ' , N.DATA 

STOP 
END XF 

CALL LSQ_POI,T2( P.COEF. TWORK, WKARRAYl, N.DATA) 

DO I = INDEXl , IW^f'^-'t^fj^lS^ 
IF {HMFLAG(I).NE.l.) THEN 

AZWORK(I) - F COEFd) 

lizWORKilT = AZWORK(I) * TIMNOM(I) - P_COEr(J) 
END DO 

^^^IzWORKd) - AZMEAS(I) 
END IF 
END DO 

IF ( I SEGMENT . EQ - "SEGMENT ) THEN 

DO I = NWEIGH*! . NAPEX+NEXTRA 

AZWORK(I) = P 

izWoiK(I) = AZWOPK(I) * TIMNOM(I) * P_COEF(J) 
END DO 
END DO 
END IF 
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END DO ! I SEGMENT 



cllllll'smoothinrranird^t^ 

""300 PRINT % 'Smooth range ? l-integ rr, 2«average' , RNGOPTION 

READ(5,*) RNGOPTION 

IF ( RNGOPTION. EQ.l) THEN 
RGWORK(J. ) * RGMEAS(.l) 

END DO 

ELSE IF ( RNGOPTION -EQ. 2) THEN 

CONTINUE 
ELSE 

GO TO 300 
END IF 



cI"II"smoot.hing data with moving *^®^*9e^«chnique 

a,N AV WEIGHT) 



PRINT 250. (AV WEIGHT5{T) UN AV WEIGHT) ^ 
PRINT ' TAR> Chance weights T (Tsyes, Osno) : 
READ ( 5 , * ) I 

IF (I.NE.ll GO TO 290 m wpthht 

PRINT ' TAR> Number of weights ? :'r N_AV_WEIGHT 

READ(5,*) I ^ 

IF ( I .LE. 0) GO TO 290 

N AV WEIGHT = I - . » 

dtjtmt * ' TAR> Enter we»-ants : 

PRIN? 25 0, TaO WEIGHTSd* ; 1 = 1, NAV WEIGHT) ^- , , . 
250 f6rHATT7X, 10{lX,F5.2),-9(7,7X,10(lX,r5.2)) 

REAI>(5,* ) AV WEIGHTS 
290 CONTINUE 

NPOINTS = NAPEX * NEXTRA 

DO I « 1 . NPOINTS 

RGWORK(I) - WKARRAiM(I) 
END DO 
END IF 
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C CALL HOV.AVERAGE( HKAKRAYl ,_RK«ORK, HPOIHTS^ 
? DO l'= 1 . N^O^lSlRM.TKt^ " 

END DO ^ 
EVD DO 



rlll ^(CINT«0-1)^ 




A2W0RK ( 



ELWORK(I) = ELWORR(A» » 
END DO 

elllll-'plorprocesserweosur^^^ 

""'do I = HAPEX*HEXTR|-1 . COUNT-1 
RGWORK(I) - SilJIti I 

ELWORK(I - f"StlfI» 
AJWORKd) - A2KEAS{I» 

END DO 



SUBSTITUTE SHEET 



wo 92/19928 ^ ^ PCr/US92/03377 



109 

CALL 5AVEDATA( RRWORK(T), DATLEN, I, 102) • 
CALL SAVEDATA ELWORK ( J. ) ♦ DTORAD , DATLEN, I, 103 
CmIl SAvi5ATA( izWORK(l)*DTORAD, DATLEN, I. 104) 
END DO 

CALL PLOTROUTINE 

^ Ask if user wanrs t-.j modify the number of segments 

IF (PDFLAG.EQ. 1 • ) THEN 
^00 ^^^^^ ^^FORMATdX, 'Modify number of time segments ? (l:yes, 0:no)') 

?f^?i'EQ^f GOTO ion 
END IF 

Q Replace real measurements by approximation 

DO I = 2 , NAPEX+NEXTPA 

PGMEAS ( I ) .» RGWORK ( 1 ) 

RRMEAS(I) = RRWORK(T ) 

ELMEASd) = ELWORKCD 

AZMEASd) = AZWORKfH 
END DO. 

RETURN 



' SUBROUTINE LIMIT_WIND( WIND, ALT. NPOINTS , MAX_SLOPE ) 

C This routine limits the variation of the wind to MAX_SLOPE. 

C For i«2,NP0INTS, the output WIND verifies: 

C abs( (WlND(i)-WIND(i-l)) / ( ALT( i) -ALT( i-1) ) ] MAX.SLOPE 

IMPLICIT NONE 
INTEGER MPOINTS 



SUBSTITUTE SHEET 



PCr/US92/03377 

WO 92/19928 



110 

REM, WIND(NPOIKTS), ALTf MPOXNTS) , JttX_SLOPE 
pSE'^IeLTA, MAXlDELTA 
n.. I = 2 , NPOINTS 

'END IF 

END DO 

^-^^'^^ 
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WHAT IS CLAIMED IS: 

1. A trajectory analysis radar system 

comprising: 

a radar antenna disposed proximate a gun barrel 
to track a projectile fired from the barrel along an 
5 actual trajectory of motion; 

a radar transmitter coupled to communicate a 
sequence of radar signal pulses to the antenna after the 
projectile is fired from the gun barrel and at least until 
the projectile reaches a zenith in the trajectory of 
10 motion; 

a radar receiver coupled to receive from the 
antenna radar signals reflected from the projectile and 
generate signals indicative of the position and velocity 

of the projectile; 

15 a signal processor coupled to receive the 

position and velocity indicative signals from the radar 
receiver and generate position and velocity data for each 
of a plxirality of different points along the trajectory of 
motion of the projectile; and 

2 0 a data processor coupled to receive the position 

and velocity information from the signal processor and 
using the received information to update atmospheric 
tables for use in a subsequent firing. 
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2 A trajectory analysis radar system 
according to claim 1 wherein the radar receiver is a 
Doppler effect type of receiver. 

3. A trajectory analysis radar system 
5 comprising: ^^^^ ^^^^^ ^^^^^^^^ ^^^^^^^^^ ^ ^ 

to traclc a projectile fired from the barrel along an 

actual trajectory of motion; . 

a radar transmitter coupled to conmunxcate a 
sequence of radar signal pulses to the antenna after the 
prl^ectile is fired from the gun barrel and at least unt.l 
the projectile reaches a zenith in the trajectory of 

motion; . ^ ^. 

a radar receiver coupled to recexve from the 

antenna radar signals reflected from the P^-^^^^^^ 
generate signals indicative of the positxon and velocity 

of the projectile; . 

a signal processor coupled to recexve the 
position and velocity ind^ative signals from the radar 
receiver and generate position and velocity data for each 
of a plurality of different points along the trajectory of 
motion of the projectile; and 

a data processor coupled to recexve the posxtxon 
and velocity information from the signal processor and 
using , the received information to update atmospherxc 
tables for use in a subsequent firing, the data processor 
assuming a current atmospheric model and repeatedly 
calculating a derived trajectory using the current 
analytical model, deriving error values representing 
differences between the actual trajectory and the derxved 
trajectory, and using the error values to correct the 
current atmospheric model to cause the current atmospheric 
model to converge toward an accurate representation of 
atmospheric conditions down range of the antenna. 

4 A trajectory analysis radar system 
according to claim 1 wherein the entire radar system is 
mounted on a weapon firing the projective. 
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5. A trajectory analysis r adar system for 
analyzing a trajectory of a projectile fired from a gun at 
a firing position, the system comprising: 

a radar system having an antenna mounted on the 
5 gun, the radar system tracking at least a portion of the 
trajectory of the projectile and generating information 
representing the tracked trajectory portion; 

a signal processing system receiving the 
trajectory representing information generated by the radar 
10 system and converting said information to numeric 
coordinate data representing the trajectory of the 

projectile; and 

a data processing system coupled to receive the 
numeric coordinate data from the signal processing system 
15 and generate in response thereto an atmospheric model 
accurately representing atmospheric conditions through 
which the trajectory of the projectile passes. 

6. A mobile weapon comprising: 

a gun firing a pro(Jectile through a trajectory; 

a radar system mounted on the weapon and 
generating radar data representing the actual trajectory 
5 of the projectile; 

a signal processing system mounted on the 
weapon, the signal processing system receiving the radar 
data and converting the radar data to coordinate based 
numeric data representing the trajectory of the 
10 projectile; and 

a fire control data processor mounted on the 
weapon and receiving the coordinate based numeric data 
from the signal processing system, the fire control data 
processor analyzing the coordinate based numeric data to 
15 generate in response thereto an atmospheric model 
representing atmospheric conditions along the trajectory 
of the projectile. 
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7. A mobile weapon comprising: 
a gun firing a projectile through a trajectory; 
a radar system mounted on the weapon and 
generating radar data representing the actual trajectory 

5 of the projectile; 

a signal processing system mounted on the 
weapon, the signal processing system receiving the radar 
data and converting the radar data to coordinate based 
numeric data representing the trajectory of the 

10 projectile; and 

a fire control data processor mounted on the 
weapon and receiving the coordinate based numeric data 
from the signal processing system, the fire control data 
processor analyzing the coordinate based numeric data to 

15 generate in response thereto an atmospheric model 
representing atmospheric conditions along the trajectory 
of the projectile, the fire control data processor 
including means for establishing a current set of 
atmospheric parameters and then repeatedly generating a 

20 derived trajectory, determining error differences between 
the actual trajectory and Ge derived trajectory, and 
using the error differences to update the current set of 
atmospheric parameters until the current set of 
atmospheric parameters converges to a model that 

25 accurately represents actual atmospheric conditions. 
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8. A mobile weapon according to claim 7 
wherein the establishing means uses a down range component 
of projectile position and velocity error to generate a 
correction value for a down range component of wind 
5 velocity, uses a cross range component of projectile 
position and velocity error to generate a correction value 
for a cross range component of wind velocity and uses an 
elevation component of position and velocity error to 
generate a correction value for air density. 



9. A mobile weapon according to claim 8 
wherein said atmospheric model represents components of 
down range and cross range wind and air density at a 
plurality of different elevation levels separated by no 
S more than 1000 foot elevation intervals. 
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10. A -thod of d.rivin, a true a»ospheri= 
u.in, » acruaX projactil. «a3.«=ry ==-pr«.n, rh. 

steps of: . 

establishing a current atmospheric modal, 
5 deriving a trajectory o£ the projectile based 

upon the current atmospheric model! .-tual 
determining error differences between th. actual 

projectile trajectory and th. derived «aj«ctory, a^^d 

correcting the curr«.t at«o.ph«r« model in 
10 response to th. determined «ror differences so a, to tend 

to reduce the error di£f«:enc.. between a derived 

trajectory based upon the curr«.t atmospheric model and 

the actual projectile trajectory. 

11. A method of deriving a true atmospheric 
model according to' claim 10 further 

rsoeating the steps otr deriving, determining and 
"^!=Sg until the current atmospheric model rapresents 
5 Ic^al .t!=.ph.ric conditions with a desire., accuracy. 

12. A method of deriving a true atmospheric 
model according to claim 11 wherein th. step, of der^ing. 
determining and corr.cting ar. repeated exactly once. 

13. X method of deriving a true atmospheric 
model according to claim 10 wherein th. ^""^^^^^ ^^^^ 
represents down range and cross rang, wmd valoc-^V 
ei^ d««ity at .l.vation intervals no greater than 

5 looo feet. 
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14. A method of firing a gun of a weapon system 

comprising the steps of: 

firing an initial round of a projectile; 

trac}cing the initial round with a radar system 
5 mounted on the weapon system to determine an actual 
trajectory of the initial round; 

determining in response to the actual trajectory 
an atmospheric model representing atmospheric conditions 
in the vicinity of the weapon system; 
10 aiming a gun of the weapon system in response to 

the determined atmospheric model; and 

firing at least one round from the gun that is 
aimed in response to the determined atmospheric model. 



15. A method of firing according to claim 14 
wherein the determining step further includes the steps 
of: 

establishing a current atmospheric model; 
5 calculating a projectile trajectory using the 

atmospheric model; 

determining error differences between the actual 
trajectory and the calculated trajectory; and 

correcting the current atmospheric model in 
10 response to the error differences to tend to cause the 
current atmospheric model to more accurately represent 
actual atmospheric conditions encountered by the initial 
round of a projectile. 
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16. A mobile weapon comprising: 
a gun firing a projectile through a trajectory; 
a radar system mounted on the weapon and 
generating radar data representing the actual trajectory 

5 of the projectile; 

a signal proneasing system mounted on the 
weapon, the signal processing system receiving the radar 
data and converting the radar data to coordinate based 
numeric data representing the trajectory of the 

10 projectile; and 

a fire control data processor mounted on the 
weapon and receiving the coordinate based numeric data 
•from the signal processing system, the fire control data 
processor analyzing the coordinate based numeric data to 

15 generate in response thereto an atmospheric model 
representing atmospheric conditions along the trajectory 
of the projectile, the fire control data processor 
establishing a current set of atmospheric parameters and 
then repeatedly generating a derived trajectory xn 

20 response to the current set of atmospheric parameters, 
determining error dilfferences between the actual 
trajectory and the derived trajectory, and using the error 
differences to update the current set of atmospheric 
parameters until the current set of atmospheric parameters 

25 converges to a model that accurately represents actual 
atmospheric conditions. 
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